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EXECUTIVE  SUMMARY 


OBJECTIVE 

The  purpose  of  this  study  was  to  upgrade  the  emergency  response  atmospheric  modeling  capabilities 
at  Vandenberg  Air  Force  Base  with  state-of- the- science  three-dimensional  forecast  models  adopted  on 
an  advanced  scientific  workstation. 

BACKGROUND 

While  employed  at  Los  Alamos  National  Laboratory,  Principal  Investigator  (PI)  of  this  study  was 
PI  of  the  project  entitled  “Three-Dimensional  Modeling  of  Rocket  Propellant  Dispersion,”  which  was 
sponsored  by  the  U.S.  Air  Force  Engineering  and  Services  Center,  Tyndall  Air  Force  Base,  FL.  The 
purpose  of  that  project  was  to  demonstrate  the  feasibility  of  using  a  desk-top  computer  to  operate  three- 
dimensional  hydrodynamic  and  diffusion  models  to  describe  the  transport  and  dispersion  of  atmospheric 
pollutants  at  VAFB.  To  accomplish  the  objective,  HOTMAC  (Higher  Order  Turbulence  Model  for 
Atmospheric  Circulation)  and  RAPTAD  (RAndom  Puff  Transport  and  Diffusion)  models  were  modified 
to  run  on  an  engineering  workstation  computer.  The  project  was  completed  in  April  1990. 

HOTMAC  and  RAPTAD  are  significantly  different  from  any  AF  models  currently  used  at  Van¬ 
denberg  Air  Force  Base  (VAFB).  HOTMAC  is  a  prognostic  model  and  solves  a  set  of  time-dependent 
physical  equations,  such  as  conservation  equations  of  momentum,  internal  energy,  mixing  ratio  of  water 
vapor,  and  turbulence  variables.  HOTMAC  forecasts  three-dimensional  distributions  of  wind  speed, 
wind  direction,  temperature,  and  moisture.  The  current  AF  models  use  a  single  wind  value  or  a  vertical 
profile  of  wind  speed  and  direction  at  a  single  location.  Clearly,  a  single  value  or  a  smgle  profile 
of  wind  cannot  provide  realistic  wind  distributions  for  accurate  transport  and  diffusion  computations, 
because  the  terrain  and  meteorology  in  the  VAFB  area  are  highly  heterogeneous.  HOTMAC  provided 
to  RAPTAD  both  mean  and  turbulence  variables  to  simulate  transport  and  diffusion  processes  of  air¬ 
borne  materials.  Only  a  few  mesoscale  atmospheric  models  can  forecast  three-dimensional  variations 
of  atmospheric  turbulence. 


The  HOTMAC/RAPTAD  model  results  were  compared  with  data  collected  during  the  Mountain  Iron 
Diffusion  experiments  (Reference  4),  which  were  conducted  at  VAFB  in  1965  and  1966.  HOTMAC 
and  RAPTAD  surface  concentration  predictions  of  fluorescent  particles  were  found  to  be  at  least  as 
good  as  those  obtained  by  diagnostic  models  where  wind  data  were  available,  and  are  potentially  far 
better  where  wind  data  are  not  available. 

Yamada  and  Bunker  (Reference  13)  showed  that  it  was  feasible  to  operate  HOTMAC  and  RAPTAD 
on  a  desk-top  computer.  HOTMAC  took  4  hours  11  minutes  and  22  hours  10  minutes  CPU  time, 
respectively,  on  a  Sun  4/1 10  workstation  and  a  Micro Vax  2000  computer  for  a  28-hour  forecast  with  21 
X  25  X  16  grid  points.  RAPTAD  took  26  minutes  and  3  hours  45  minutes  CPU  time,  respectively,  on 
a  Sun  4/110  workstation  and  a  Micro  Vax  2000  computer  for  a  20  hour  simulation.  Significantly  faster 
(approximately  7  times  faster  than  a  Sun  4/110)  workstations  have  become  available  since  the  previous 
project  was  completed.  The  CPU  times  quoted  above  become  much  shorter  with  these  new  workstations. 

The  increased  affordability  and  portability  of  workstations  have  opened  the  door  to  upgrading 
toxic -hazard  modeling  capabilities  for  emergency  response  management.  The  HOTMAC  and  RAPTAD 
modeling  system  will  be  a  useful  enhancement  to  current  launch  operation  planning  and  emergency 
response  management  capabilities  at  VAFB. 

TASKS 

The  following  tasks  were  performed  and  are  described  in  detail  in  Section  III. 

Task  1.  Develop  a  method  to  incorporate  National  Meteorological  Center’s  (NMC)  Limited- Area 
Fine  Mesh  (LFM)  and/or  Nested  Grid  Model  (NGM)  weather  forecast  into  HOTMAC. 

Task  2.  Develop  a  method  to  restart  HOTMAC  prediction,  and  maintain  on  disk  predictions  for 
the  next  24  hours. 

Task  3.  Develop  a  method  to  combine  upper-air  wind-sounding  and  tower  data  taken  at  VAFB  with 
HOTMAC  forecast  winds,  so  that  RAPTAD  can  use  the  composite  (forecast  and  observed)  winds. 

Task  4.  Add  to  the  model  the  physics  necessary  to  forecast  the  formation  and  dissipation  of  fog. 


VI 


Task  5.  Develop  an  innovative  method  to  treat  plumes  that  are  not  neutrally  buoyant. 

Task  6.  Develop  a  method  to  predict  concentration  variances,  which  are  used  to  determine  how  a 

required  confidence  level  affects  the  extent  of  a  critical  area. 

Task  7.  Perform  model  verification  runs  with  a  nested  grid  version  of  HOTMAC  on  a  workstation. 

Task  8.  Upgrade  computer  capabilities  used  for  emergency  response  modeling  at  VAFB. 
METHODOLOGY 

Section  III.A.l  reviewed  briefly  the  model  equations,  boundary  conditions,  and  numerical  schemes 
used  in  HOTMAC  and  RAPTAD.  HOTMAC  is  a  mesoscale  atmospheric  model  that  can  forecast  three- 
dimensional  distributions  of  wind  speed,  wind  direction,  turbulence,  temperature,  and  water  vapor.  The 
basic  equations  for  HOTMAC  are  the  conservation  equations  for  mass,  momentum,  internal  energy, 
mixing  ratio  of  water  vapor,  and  turbulence  kinetic  energy.  Model  physics  in  HOTMAC  were  improved 
considerably  in  order  to  simulate  evaluation  of  cloud  formation  and  precipitation  processes.  Interaction 
among  radiation  transfer,  cloud  development,  and  turbulent  motion  was  included.  Time  variations  of 
large-scale  wind  speed,  wind  direction,  potential  temperature,  and  mixing  ratio  of  water  vapor  can  be 
incorporated  into  HOTMAC  by  a  four-dimensional  data  assimilation  method. 

RAPTAD  is  a  Lagrangian  puff  code  which  uses  the  Monte  Carlo  statistical  diffusion  process.  The 
center  location  and  standard  deviations  of  concentration  distribution  for  each  puff  are  computed  by  use 
of  wind  and  turbulence  which  are  modeled  by  HOTMAC.  Concentration  at  any  location  is  computed 
by  the  summation  of  the  concentrations  contributed  by  all  of  the  puffs.  RAPTAD  can  be  used  under 
extreme  conditions  with  highly  heterogeneous  wind  and  turbulence  distributions  where  a  conventional 
Gaussian  plume  model  may  fail.  Tower  data  may  be  incorporated  into  RAPTAD  in  addition  to  the 
wind  distributions  provided  by  HOTMAC.  Buoyancy  effects  of  plumes  (positive  and  negative)  were 
also  included  in  RAPTAD. 

RESULTS 

Task  1:  Four-Dimensional  Data  Assimilation 
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To  test  the  simulation  capabilities  of  HOTMAC,  we  selected  data  from  Phase  I  WIND  data 
(Reference  18).  Horizontal  wind  components  were  nudged  to  the  target  wind  components  computed 
from  observed  winds,  the  Coriolis  parameter,  and  a  nudging  parameter  used  in  the  equations  of  motion. 

The  Phase  I  was  charactenzed  by  typical  summertime  anticyclomc  conditions  over  the  Sacramento 
River  Valley  in  northern  California.  Southerly  flow  (upslope  wind)  during  the  day  and  northerly  flow 
(downslope  wind)  at  night  were  the  characteristic  flow  patterns. 

Comparisons  between  observations  and  simulations  of  the  wind  directions  and  speed,  temperature, 
and  dew  point  showed  good  agreement,  both  spatially  and  temporally,  throughout  entire  24-hour  periods. 
Thus,  a  single  station  sounding  could  be  used  to  represent  large-scale  meteorological  conditions  over 
the  model  domain. 

Statistical  parameters,  including  the  mean  and  standard  deviation  of  both  the  observed  and  simulated 
variables,  the  systematic  and  unsystematic  components  of  the  root-mean-square  differences  as  well  as 
the  total  root-mean-square  difference  (RMSDs,  RMSDU,  and  RMSD,  respectively),  and  the  agreement 
measures,  were  calculated  hourly  using  the  surface  station  data.  These  statistical  parameters  were 
calculated  for  three  different  wind  fields  simulated  by  nudging  to  the  target  winds,  with  no  nudging, 
and  persistence.  The  wind  fields  obtained  by  nudging  to  target  winds  yielded  better  agreement  in  wind 
direction  and  speed  when  compared  to  the  two  other  wind  fields. 

Wind  Distributions  at  6  Meters  agl  at  1500  LT  on  Day  233.  Similar  to  Figures  37  and  43,  But 
Without  a  Nested  Grid 

Task  2:  Continuous  Forecast 

A  complex  scheme  was  developed  to  operate  HOTMAC  and  RAPTAD  interactively.  HOTMAC 
runs  continuously  and  always  maintains  24— hour  forecast  data.  The  interactive  scheme  dictates  that  as 
soon  as  RAPTAD  starts  computation,  HOTMAC  halts  its  forecasting. 

Once  a  24-hour  forecast  is  completed  and  the  forecast  data  are  written  to  a  disk,  HOTMAC  "sleeps" 
until  the  next  computation  time.  HOTMAC  "wakes  up"  at  the  clock  hour  and  begins  to  advance  the 
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forecast  one  hour. 


In  addition  to  using  the  default  mode’s  24-hour  forecast,  the  user  can  extend  the  range  of  the 
forecast  period  manually,  up  to  48  hours  or  longer,  for  planning  operations. 

Task  3:  Tower  Data  Input 

Winds  measured  by  upper-air  soundings  and  at  towers  were  combined  with  the  winds  forecasted  by 
HOTMAC.  A  simple  Hr2  weighting  method  was  used  to  compute  the  wind  at  a  puff  center.  The  winds 
at  puff  centers  (obtained  from  tower  winds)  were  then  combined  with  winds  computed  by  HOTMAC. 
Tower  winds  are  weighted  based  on  the  distance  between  a  puff  center  and  the  nearest  tower  site.  If 
a  puff  is  located  at  a  tower  site,  then  the  tower  wind  is  used  for  puff  velocity.  The  weight  function 
decreases  exponentially  as  the  distance  between  a  puff  and  a  tower  increases.  Such  an  assumption  is 
reasonable,  particularly  for  airflows  over  complex  terrain,  where  wind  distributions  are  largely  localized. 

Task  4:  Simulation  of  Fog 

Fog  and  low  stratus  clouds  are  frequently  observed  at  VAFB,  and  affect  the  heat  energy  balance  at 
the  ground.  The  interaction  of  fog,  solar  heating,  long-wave  radiation  heating/cooling,  and  turbulence 
mixing  is  complex  and  is  not  yet  well  understood.  A  simulation  of  fog  consists  of  modeling  condensation, 
radiation  transfer,  and  precipitation  microphysics.  Condensation  processes  were  simulated  by  an 
ensemble  cloud  model  which  is  described  in  detail  in  Section  ni.D.l. 

The  terms  for  the  rate  of  condensation  are  purposely  eliminated  by  introducing  the  liquid  water 
potential  temperature  and  mixing  ratio  of  total  water.  In  order  to  recover  the  potential  (or  absolute) 
temperature  and  the  mixing  ratios  of  water  vapor  and  liquid  water,  Gaussian  cloud  relations,  proposed 
by  Sommeria  and  Deardorff  (Reference  33)  and  Mellor  (Reference  34),  are  used.  The  present  method 
has  been  applied  to  simulations  of  the  BOMEX  data  (Reference  35),  and  GATE  data  (Reference  36). 

Solar  and  long-wave  radiation  play  important  roles  in  the  formation  and  dissipation  of  clouds  and  fog. 
Hanson  and  Derr  (Reference  39)  proposed  a  parameterized  solar  radiation  scheme  whose  parameters 
were  obtained  by  curve-fitting  to  numerical  radiative-transfer  results  using  the  ATRAD  narrow-band 
model  (Reference  40).  This  parameterized  method  is  simple,  yet  reproduced  solar  flux  profiles  within  a 
single  cloud  layer  which  were  in  good  agreement  with  the  numerical  results  obtained  from  ATRAD. 


IX 


Motivated  by  the  success  of  the  solar  radiation  scheme,  Hanson  and  Derr  proposed  an  infrared 
(IR)  radiation  scheme  expressed  by  exponential  functions  whose  decay  parameters  were  determined  by 
the  emissivity  methods  (Reference  41).  This  method  reproduced  IR  flux  profiles  within  layered  clouds 
which  were  in  good  agreement  with  the  numerical  results  and  observations  for  thick  clouds  (-800  meters). 
However,  the  parameterization  overestimated  the  flux  decrease  at  the  cloud  top  and  underestimated  the 
cloud-base  warming  compared  to  the  numerical  model  for  thin  clouds  (-300  meters). 

We  adopted  the  Hanson  and  Derr  parameterization  scheme  because  of  its  simplicity  and  because  it 
can  produce  flux  profiles  which  are  in  good  agreement  with  observations  and  numerical  model  results. 

We  adopted  the  precipitation  microphysics  proposed  by  Nickerson  et  al.  (Reference  43).  This  hy¬ 
drological  model  is  substantially  more  complex  than  Kessler’s  (Reference  44)  classical  parameterization 
for  the  treatment  of  microphysics.  The  model  is  based  on  a  more  realistic  log-normal  distribution  than 
the  Marshall-Palmer  distribution  in  Kessler’s  parameterization.  The  mean  diameter  of  raindrops  in  a 
given  volume  is  estimated  from  equations  for  the  mean  rain  water  mixing  ratio  and  the  raindrop  number 
concentration.  The  terminal  velocity  of  raindrops  in  each  volume  is  calculated  as  a  function  of  the 
mean  raindrop  diameter,  thus  the  terminal  velocities  in  the  model  are  expected  to  be  more  accurate  than 
those  obtained  by  a  simple  parameterization. 

Maritime  stratus  clouds  were  simulated  by  using  a  one-dimensional  version  of  HOTMAC  (Reference 
20).  Integration  initiated  at  0000  LT  (Local  Tune)  on  day  200  (July  19)  and  continued  for  48  hours.  The 
results  for  the  second  24-hour  simulations  are  presented  here.  Clouds  were  thick  during  the  nocturnal 
period  and  occupied  the  layer  between  1200  meters  and  1600  meters  above  the  sea  surface.  As  the 
short-wave  solar  radiation  heating  increased,  clouds  dissipated  from  the  lower  part  and  clouds  became 
as  thin  as  100  meters  by  1800  LT.  As  the  solar  heating  subsided,  the  cloud  thickness  increased  again. 

A  one-dimensional  version  of  HOTMAC  was  used  to  simulate  time  evolution  of  fog  over  a 
horizontally  homogeneous  terrain.  The  results  were  compared  with  data  taken  at  Cabaw  meteorological 
tower  in  the  Netherlands  (Reference  31).  Fog  began  to  form  almost  immediately,  as  the  air  temperature 
decreased  due  to  cooling  at  the  ground.  The  height  of  the  fog  increased  as  the  air  temperature  near  the  fog 
top  decreased  due  to  long-wave  radiation  cooling.  The  modeled  and  observed  temperature  and  mixing 


ratio  of  water  vapor  at  1.1  meters  above  ground  level  (agl)  were  in  good  agreement  with  observations. 

Task  5:  Non- neutrally  Buoyant  Plumes 

Highly  buoyant  plumes  modify  the  wind  and  turbulence  distributions  of  the  ambient  flow.  It  is 
almost  impossible  to  parameterize  or  express  such  modifications  without  deploying  a  dynamic  plume 
model.  A  physically  correct  way  to  treat  non-neutrally  buoyant  plumes  is  to  incorporate  plume  dynamics 
into  HOTMAC.  However,  the  dynamic  plume  model  requires  considerable  computer  time  and  is  not 
practical  for  emergency  modeling.  This  is  a  dilemma,  for  one  must  choose  between  being  accurate 
(but  impractical)  and  being  practical  (but  inaccurate).  It  is  therefore  necessary  to  adopt  a  plume 
parameterization  scheme  (instead  of  a  dynamic  plume  model)  to  meet  the  time  constraints  of  an 
emergency  response  modeling  system. 

Following  Van  Dop  (Reference  45),  the  vertical  velocity  of  a  buoyant  plume  is  computed  from  the 
Langevin  equation  of  motion  for  a  homogeneous  and  stationary  turbulent  flow.  The  temperature  of  a 
buoyant  plume  is  also  assumed  to  be  computed  by  the  Langevin  equation. 

Test  simulations  were  conducted  under  neutral  conditions  (N2  =  0)  and  the  results  were  compared 
with  analytical  solutions.  Simulations  and  analytical  solutions  were  in  good  agreement.  Simulations 
were  repeated  under  stable  (N2  >  0)  and  unstable  (N2  <  0)  atmospheric  conditions.  Finally,  simulations 
were  conducted  for  a  plume  whose  initial  density  is  heavier  than  the  ambient  air. 

All  simulations  were  conducted  with  an  integration  time  step  of  0.1  seconds.  Accuracy  decreased 
rapidly  with  increases  in  the  integration  time  step.  RAPTAD  uses  a  time  step  of  10  seconds,  but  a  small 
time  step  was  necessary  to  simulate  initial  plume  rise  (drop)  accurately.  Therefore,  a  time  step  which 
increased  linearly  with  time  was  used:  time  step  was  0.1  seconds  initially  and  became  10  seconds  at 
100  seconds  after  the  release. 

Task  6:  Concentration  Variance 

Relatively  short  time- averaging  values  are  required  for  predicting  concentrations  of  toxic  materials. 
Such  values  normally  exhibit  great  variations  in  time  and  space.  Thus,  predictions  of  not  only  the 
mean  but  also  the  variance  of  concentrations  are  essential  to  determine  the  extent  of  a  critical  area 
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where  concentration  values  exceed  a  given  limit.  The  size  of  a  critical  area  becomes  larger  if  a  higher 
confidence  level  is  required. 

Only  a  limited  number  of  concentration  fluctuation  measurements  are  available  for  comparison. 
The  best  measurements  were  obtained  by  using  wind  tunnels.  Preliminary  comparisons  indicate  that 
our  simulations  are  qualitatively  in  good  agreement  with  wind  tunnel  data.  Both  standard  deviations, 
and  ratios  of  standard  deviation  to  mean,  showed  characteristics  similar  to  those  found  in  measurements 
where  measured  mean  wind  and  turbulence  were  used  in  RAPTAD  computation.  No  effort  was  made 
to  compare  simulations  with  atmospheric  data  because  we  are  not  aware  of  such  data. 

Task  7:  Nested  Grid 

A  nested  grid  version  of  HOTMAC  had  been  used  with  a  Cray  supercomputer,  but  it  had  never  been 
used  with  a  workstation  because  it  would  require  too  much  computer  time.  However,  with  the  rapid 
advancement  of  workstation  performance,  we  are  now  able  to  run  a  nested  grid  version  of  HOTMAC 
on  a  workstation.  Test  simulations  have  been  performed  to  demonstrate  systematically  how  nested  grids 
improve  the  simulation  of  wind  distributions.  A  two-way  nesting  method  was  used,  the  outer  grid 
provided  boundary  values  to  the  inner  grid,  and  the  inner  grid  updated  outer  grid  values  with  new 
values  at  grid  points  common  to  both  grids. 

Simulations  were  conducted  with  and  without  nested  grids.  Three  nested  grids  were  used  in  a 
control  run.  Wind  distributions  at  1500  LT  for  a  single  grid  case  were  similar  to  those  for  the  control 
run  because  strong  mixing  in  the  vertical  direction,  due  to  turbulence,  resulted  in  almost  uniform  wind 
distributions.  On  the  other  hand,  wind  distributions  at  0200  LT  were  considerably  different  from  those 
for  the  control  run,  particularly  in  the  area  where  Grid  3  was  nested. 

Task  8:  Computer  Capability  Upgrade 

Based  on  discussions  between  the  personnel  at  VAFB  and  YSA,  the  following  computer  hardware 
and  software  were  selected: 

1.  A  Sun  SPARCstation  10  with  32  Mb  RAM,  19"  color  monitor,  400  Mb  internal  SCSI  disk,  1.44 
Mb  3t"  internal  floppy  disk. 
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2.  A  1.2  Gb  external  disk. 

3.  A  150  Mb  tape  drive. 

4.  A  CD  drive. 

5.  Sun  FORTRAN  compiler. 

6.  Sun  C  compiler. 

7.  NCAR  Graphics. 

Hardware  and  software  were  installed  and  tested  at  YSA,  and  delivered  to  VAFB. 

CONCLUSION 

•  It  is  feasible  to  operate  on  a  workstation,  a  nested-grid,  three-dimensional  atmospheric  model, 
HOTMAC,  and  a  three-dimensional  dispersion  model,  RAPTAD,  to  forecast  the  transport  and 
diffusion  of  airborne  materials  at  VAFB. 

•  The  computer  capabilities  for  emergency  response  applications  at  VAFB  were  upgraded.  HOT¬ 
MAC  and  RAPTAD  were  installed  and  tested  on  a  Sun  SPARCstationlO. 

•  HOTMAC  was  modified  to  run  continuously  and  store  wind  and  turbulence  predictions  and  data 
for  the  next  24  hours.  When  an  emergency  occurs,  RAPTAD  can  be  used  immediately  with  the 
wind  data  stored  on  a  disk.  The  RAPTAD  computation  is  much  faster  than  that  of  HOTMAC. 
This  approach  meets  better  the  time  constraints  of  emergency  situations. 

•  A  method  was  developed  to  integrate  large-scale  weather  data  into  HOTMAC.  The  weather  data 
are  used  to  initialize  and  correct  HOTMAC  forecasts.  A  four-dimensional  data  assimilation 
method  was  used. 

•  A  method  was  developed  to  predict  concentration  variances  which  can  be  used  to  estimate 
uncertainties  associated  with  predictions. 

•  Model  physics  necessary  to  simulate  the  evolution  of  fog  formation  and  dissipation  processes 
were  added.  Fog  is  frequently  observed  at  VAFB  and  affects  the  heat  energy  balance  at  ground 
level. 

•  Positive  and  negative  buoyancy  effects  of  plumes  were  incorporated  to  RAPTAD. 
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recommendations 


•  Become  familiar  with  background  theories  and  operating  procedures  of  HOTMAC  and  RAPTAD. 
Increase  experience  by  running  models  under  different  weather  conditions. 

o  Consider  customizing  pre-  and  post-  processors  of  HOTMAC  and  RAPTAD  to  meet  specific 
needs  at  VAFB.  The  current  versions  of  HOTMAC  and  RAPTAD  are  designed  for  general 
applications. 

3  Continue  to  upgrade  computer  capabilities  at  VAFB.  Computer  hardware  technology  is  expected 
to  advance  further  and  to  allow  much  more  sophisticated  models  to  be  operated  in  much  less 
time  than  previously  considered  possible. 

®  Continue  to  improve  model  physics  in  HOTMAC  and  RAPTAD.  These  codes  have  model 
physics  which  are  considered  to  be  the  state-of-the-science  but  still  are  far  from  fully  replicating 
the  complex  physics  of  the  real  atmosphere.  Particularly,  we  recommend  strongly  that' the 
new  model  physics  added  here  (precipitation  microphysics,  radiation  transfer  in  fog  and  clouds, 
and  positive  and  negative  plume  buoyancy  effects)  be  tested,  by  comparing  simulations  with 
observations. 
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SECTION  I 


INTRODUCTION 


A.  OBJECTIVE 

The  objective  of  this  effort  was  to  upgrade  the  emergency  response  atmospheric  modeling 
capabilities  at  VAFB  with  state-of-the-science  three-dimensional  forecast  models  adopted  on  an 
advanced  scientific  workstation.  Air  Force  (AF)  models  designed  to  calculate  the  atmospheric 
transport  and  dispersion  of  pollutants  at  Vandenberg  Air  Force  Base  (VAFB)  include  the  Rocket 
Exhaust  Effluent  Dispersion  Model  (REEDM)  system  (Reference  1;  Reference  2)  and  the  Ocean 
Breeze/Dry  Gulch  (OB/DG)  equation  (Reference  3),  Mountain  Iron  (MI)  equation  (Reference  4; 
Reference  5),  and  Sudden  Ranch  (SR)  equation  (Reference  6).  REEDM  was  designed  to  model 
the  transport  and  dispersion  of  buoyant  plumes  of  rocket  exhaust  during  normal  launch,  and 
fireball  rise  and  subsequent  dispersion  in  the  event  of  an  explosive  accident.  OB/DG,  MI,  and 
SR  are  empirical  equations  designed  to  calculate  toxic-hazard  corridors  resulting  from  accidental 
spills  of  toxic  chemicals.  These  models  assume  that  the  wind  field  over  the  geographical  area 
of  interest  can  be  adequately  represented  by  the  wind  speed  and  direction  measured  at  a  single 
point  (or,  in  the  case  of  REEDM,  a  vertical  profile  of  wind  speed  and  direction  over  a  single 
site).  These  models  also  assume  that  the  toxic  plume  or  puff  released  to  the  atmosphere  is 
uniformly  distributed  around  a  center  line  determined  from  the  wind  direction  measured  at  the 
single  observation  site. 

While  the  model  and  equations  listed  above  have  been  acceptable  in  the  past,  the  complex 
terrain  and  meteorology  at  VAFB,  increased  amounts  of  toxic  chemicals  stored  there,  increases 
in  the  surrounding  civilian  population,  higher  anticipated  launch  rates,  and  the  reduction  in 
personal  exposure  limits  for  the  AF  chemicals  of  interest  (nitrogen  tetroxide  and  Aerozine-50) 
require  more  accurate  methods.  Current  efforts  to  improve  VAFB’s  modeling  capabilities  include 
development  of  the  Air  Force  Toxic  (AFTOX)  Chemical  Dispersion  model  (Reference  7)  by  the 
Air  Force  Geophysics  Laboratory  (AFGL).  AFTOX,  in  its  current  form,  however,  also  assumes 
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that  the  three-dimensional  wind  field  can  be  represented  by  a  single  wind  value.  Efforts  have  been 
made  to  couple  AFTOX  with  a  surface-layer  wind-flow  model;  however,  the  complexity  of  the 
terrain  and  meteorology  at  VAFB  require  a  model  that  takes  into  account  the  three-dimensional 
wind  field  over  the  entire  VAFB  region. 

A  considerable  number  of  wind-field  models  exist  that  could  produce  three-dimensional  wind 
distributions  over  complex  terrain.  These  wind  models  are  of  two  types. 

1.  Diagnostic  wind  models  that  construct  the  wind  distribution  by  using  existing  measure¬ 
ments  plus  some  physical  constraints  of  varying  extent. 

2.  Prognostic  wind  models  that  use  a  set  of  hydrodynamic  equations  to  predict  the  wind 
and  turbulence  fields. 

Diagnostic  wind  models  are  based  on  measurements  and  interpolation.  Most  of  these 
models  adjust  initially  interpolated  wind  fields  to  satisfy  a  mass-conservation  equation.  The 
accuracy  of  the  model  results  is  greatly  dependent  upon  the  quality  and  representativeness  of 
the  measurements. 

Diagnostic  wind  models  have  been  used  almost  exclusively  for  emergency  response  manage¬ 
ment  because  they  are  simple  and  fast.  The  critical  weakness  of  this  class  of  models,  however, 
is  their  lack  of  predictive  capability.  Plume  models  based  upon  diagnostic  wind  models  must 
assume  that  the  current  wind  field  will  persist  until  it  is  updated  by  future  measurements.  The 
assumption  of  persistency  is  adequate  if  the  measurements  are  made  frequently  and  if  no  abrupt 
change  in  the  wind  field  occurs.  However,  such  conditions  limit  the  area  of  concern  to  the 
immediate  vicinity  of  the  source,  and  exclude  predictions  during  the  transitional  periods  between 
daytime  and  nighttime  conditions.  Wind  speeds,  wind  directions,  and  turbulence  change  rapidly 
during  the  transitional  periods. 

The  question  of  the  representativeness  of  measurements  further  limits  the  applications  of 
diagnostic  wind  models  to  simple  terrains,  because  complex  surface  boundaries  require  a 
prohibitively  large  number  of  wind  measurements. 
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Prognostic  models,  on  the  other  hand,  solve  a  set  of  time-dependent  physical  equations,  such 
as  conservation  equations  of  momentum,  internal  energy,  and  the  mixing  ratio  of  water  vapor. 
Prognostic  models  forecast  three-dimensional  wind  field  distributions  that  become  the  input  to 
transport  and  diffusion  models.  For  this  reason,  the  terms  “prognostic”  and  “forecast”  are  used 
interchangeably  in  this  proposal. 

At  present,  prognostic  models  have  not  been  used  in  emergency  response  systems,  because 
a  large  amount  of  computation  time  was  required.  This  situation  should  improve  considerably 
in  the  near  future,  because  the  capabilities  of  engineering  workstations  are  expected  to  advance 
at  an  astonishing  rate.  Previous  AF  efforts  to  improve  toxic-hazard  models  have  been  limited 
by  available  computer  power.  Recent  developments  in  computer  hardware  technology,  however, 
have  the  potential  to  allow  much  more  sophisticated  models  to  be  operated  in  much  less  time 
than  previously  considered  possible. 

B.  BACKGROUND 

There  have  been  considerable  efforts  to  develop  three-dimensional  wind  models  that  can 
be  used  to  simulate  wind  distributions  over  the  area  surrounding  Vandenberg  Air  Force  Base. 
A  comprehensive  review  of  the  present  emergency  modeling  system  (REEDM)  at  VAFB, 
discussions  on  candidate  wind-field  models  (both  diagnostic  and  prognostic)  for  replacement 
of  the  current  REEDM  wind-field  model,  and  recommendations  for  short-  and  long-term  future 
improvements  are  given  in  recent  reports  by  Conley  (Reference  8). 

For  example.  Figure  1  shows  the  modeled  horizontal  wind  vectors  at  4  meters  above-ground- 
level  (agl)  by  using  the  WOCSS  (Winds  On  Critical  Stream  Surfaces)  diagnostic  wind  model 
(Reference  9).  WOCSS  is  the  enhanced  version  of  SRI  Complex  developed  by  SRI  International 
(Reference  10).  An  array  of  50  X  80  X  6  with  a  500-meter  grid  spacing  and  the  terrain  data  of  a 
500-meter  resolution  were  used.  Wind  vectors  at  every  fourth  grid  point  are  shown  in  Figure  1. 

Figure  2  shows  the  modeled  wind  vectors  at  10  meters  agl  at  1310  Local  Time  (LT)  by  using 
LINCOM  (Reference  11),  which  is  a  linear,  diagnostic,  spectral,  and  potential-flow  wind  model. 
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Figure  1 .  Modeled  Horizontal  Wind  Vectors  for  the  MI87  Case  at  4  Meters  agl  from  the  WOCSS 
Diagnostic  Wind  Model.  Meteorological  Measurement  Stations  are  identified  by 
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Figure  2.  Modeled  Horizontal  Wind  Vectors  for  the  MI87  Case  at  10  Meters  agl  from  the 
LINCOM  Diagnostic  Wind  Model. 


5 


Both  WOCSS  and  LINCOM  used  the  same  tower  data,  and  wind  fields  were  adjusted  to  satisfy 
mass  consistency.  Indeed,  wind  distributions  shown  in  Figures  1  and  2  are  in  close  agreement, 
although  LINCOM  (Figure  2)  produced  slightly  more  spatial  variations  than  WOCSS  (Figure  1). 

As  mentioned  earlier,  a  diagnostic  wind  model  is  no  more  than  a  scheme  for  simple 
interpolation/extrapolation  of  observed  winds  at  specified  grid  points.  For  this  reason,  the  spatial 
variations  of  wind  vectors  in  Figures  1  and  2  are  confined  mainly  to  the  regions  where  data  are 
available.  Mass  consistency  and  other  physical  constraints  used  in  diagnostic  wind  models  affect 
the  wind  vectors  away  from  the  measurement  sites  very  little,  because  initial  wind  distributions 
are  normally  determined  by  assuming  that  the  influence  of  measurements  will  drop,  inversely 
proportionally  to  the  square  of  the  distance  from  the  measurement  sites. 

A  prognostic  wind  model,  on  the  other  hand,  is  based  on  a  set  of  physical  equations  and  does 
not  require  the  measurement  of  surface  winds  for  forecast  of  wind  and  turbulence  distributions. 
Observations  are  used  to  verify  the  predictions.  Recently,  Yamada  (Reference  12)  used  a  three- 
dimensional  prognostic  wind  model  HOTMAC,  Higher  Order  Turbulence  Model  for  Atmospheric 
Circulation,  to  simulate  diurnal  variations  of  wind  and  turbulence  distributions  in  the  Vandenberg 
area. 

HOTMAC  model  results  (Figure  3)  generally  agree  with  the  diagnostic  model  results  (Figures 
1  and  2)  in  the  area  where  wind  measurements  are  available.  Wind  vectors  in  the  prognostic 
model  are  determined  from  the  balance  between  the  pressure  force,  turbulent  mixing,  and 
Coriolis  force.  Therefore,  prognostic  models  can  provide  wind  distributions  in  regions  where 
no  measurements  are  available.  HOTMAC  vectors  (Figure  3)  are  quite  different  from  the  wind 
fields  computed  by  WOCSS  (Figure  1)  and  LINCOM  (Figure  2)  in  the  areas  where  topographic 
features  are  prominent  and  no  wind  measurements  are  available.  The  diagnostic  models  do  not 
incoiporate  thermal  pressure  forces  generated  by  the  differential  heating  over  sloped  surfaces. 
Therefore,  the  variations  of  wind  distributions  in  a  diagnostic  wind  model  must  be  provided  by 
the  wind  data  that  reflect  wind  perturbations  caused  by  surface  heterogeneity. 
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Figure  3.  Modeled  Horizontal  Wind  Vectors  for  the  MI87  Case  at  10  Meters  agl  from  the 
HOTMAC  Prognostic  Wind  Model 


C.  SCOPE 


Section  II  discusses  the  previous  study  which  is  directly  related  to  the  present  study.  Section 
III  discusses  in  detail  the  work  performed  for  each  task.  Large-scale  weather  variations  were 
incorporated  into  HOTMAC  by  using  a  technique  known  as  four-dimensional  data  assimilation. 
Model  performance  was  evaluated  by  using  several  statistical  measures,  such  as  mean  differences, 
standard  deviations,  root-mean-square  differences,  and  agreement  measures.  Section  HI  also 
discusses  newly  added  model  physics,  including  radiation  transfer  in  fog  and  clouds,  precipitation 
microphysics,  and  positive  and  negative  buoyant  plumes.  Improvements  obtainable  through  the 
use  of  nested  grids  are  also  demonstrated.  Conclusions  and  recommendations  for  future  work 
are  given  in  Section  IV. 
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SECTION  II 


PREVIOUS  STUDY 


While  employed  at  Los  Alamos  National  Laboratory,  the  Principal  Investigator  (PI)  of  this 
proposal  was  PI  of  the  project  entitled  “Three-Dimensional  Modeling  of  Rocket  Propellant  Dis¬ 
persion,”  which  was  sponsored  by  the  U.S.  Air  Force  Engineering  and  Services  Center,  Tyndall 
Air  Force  Base,  FL.  The  purpose  of  that  project  was  to  demonstrate  the  feasibility  of  using  a 
desk-top  computer  to  operate  three-dimensional  hydrodynamic  and  diffusion  models  to  describe 
the  transport  and  dispersion  of  atmospheric  pollutants  at  VAFB.  To  accomplish  the  objective, 
FIOTMAC  (Higher  Order  Turbulence  Model  for  Atmospheric  Circulation)  and  RAPTAD  (RAn- 
dom  Puff  Transport  And  Diffusion)  models  were  modified  to  run  on  an  engineering  workstation 
computer.  The  project  was  completed  in  April  1990. 

HOTMAC  and  RAPTAD  are  significantly  different  from  any  AF  models  mentioned  in  Section 
I.  HOTMAC  is  a  prognostic  model  which  solves  a  set  of  time-dependent  physical  equations, 
such  as  conservation  equations  of  momentum,  internal  energy,  mixing  ratio  of  water  vapor, 
and  turbulence  variables.  HOTMAC  forecasts  three-dimensional  distributions  of  wind  speed, 
wind  direction,  temperature,  and  moisture.  The  current  AF  models  use  a  single  wind  value  or 
a  vertical  profile  of  wind  speed  and  direction  at  a  single  location.  Clearly,  a  smgle  value  or 
a  single  vertical  profile  of  wind  cannot  yield  realistic  wind  distributions  for  accurate  transport 
and  diffusion  computations,  because  the  terrain  and  meteorology  in  the  VAFB  area  are  highly 
heterogeneous.  HOTMAC  provided  to  RAPTAD  both  mean  and  turbulence  variables  to  simulate 
transport  and  diffusion  processes  of  airborne  materials.  Only  a  few  mesoscale  atmospheric 
models  can  forecast  three-dimensional  variations  of  atmospheric  turbulence. 

The  HOTMAC/RAPTAD  model  results  were  compared  with  data  collected  during  the  Moun¬ 
tain  Iron  Diffusion  experiments  (Reference  4),  which  were  conducted  at  VAFB  in  1965  and  1966. 
HOTMAC  and  RAPTAD  surface  concentration  predictions  of  fluorescent  particles  were  found 
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to  be  at  least  as  good  as  those  obtained  by  diagnostic  models  where  wind  data  were  available, 
and  are  potentially  far  better  where  wind  data  are  not  available. 

Yamada  and  Bunker  (Reference  13)  showed  that  it  was  feasible  to  operate  HOTMAC  and 
RAPTAD  on  a  desk-top  computer.  HOTMAC  took  4  hours  1 1  minutes  and  22  hours  10  minutes 
CPU  time,  respectively,  on  a  Sun  4/110  workstation  and  a  MicroVax  2000  computer  for  a 
28-hour  forecast  with  21  X  25  X  16  grid  points.  RAPTAD  took  26  minutes  and  3  hours  45 
minutes  CPU  time,  respectively,  on  a  Sun  4/110  workstation  and  a  MicroVax  2000  computer 
for  a  20-hour  simulation.  Significantly  faster  (approximately  15  times  faster  than  a  Sun  4/1 10) 
workstations  have  become  available  since  the  previous  project  was  completed.  The  CPU  times 
quoted  above  become  much  shorter  with  these  new  workstations. 

The  increased  affordability  and  portability  of  workstations  have  opened  the  door  to  upgrad¬ 
ing  toxic-hazard  modeling  capabilities  for  emergency  response  management.  The  HOTMAC 
and  RAPTAD  modeling  system  will  be  a  useful  enhancement  to  current  emergency  response 
management  capabilities  at  VAFB. 

To  enhance  the  usefulness  of  HOTMAC  and  RAPTAD  for  operational  use  at  VAFB,  Yamada 
and  Bunker  (Reference  13)  recommended  the  following: 

1 .  Consider  upgrading  computer  capabilities  for  emergency  response  applications  at  VAFB . 
MicroVax  computers  may  not  be  fast  enough  to  operate  HOTMAC  and  RAPTAD 

usefully. 

2.  Modify  HOTMAC  to  run  continuously  and  store  wind  and  turbulence  data  for  the  next 
24  hours.  When  an  emergency  occurs,  RAPTAD  can  be  used  immediately  with  the  wind 
data  stored  on  disk.  The  RAPTAD  computation  is  much  faster  than  that  of  HOTMAC, 
and  is  thus  more  useful  in  emergency  situations. 

3.  Develop  a  method  to  integrate  tower  data  into  HOTMAC.  The  wind  data  can  be  used  to 
initialize  and  correct  the  wind  distribution  in  HOTMAC.  This  may  be  accomplished  by 
using  a  dynamic  initialization  technique  and  a  four-dimensional  data  assimilation  method. 
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4.  Develop  a  method  to  predict  concentration  variances  which  can  be  used  to  estimate 
uncertainties  associated  with  predictions.  A  short  time-averaging  value  is  required 
for  predicting  concentration  of  toxic  materials.  Such  a  value  normally  exhibits  great 
variations  in  time  and  space. 

5.  Add  model  physics  necessary  to  simulate  the  evolution  of  fog  formation  and  dissipation 
processes.  Fog  is  frequently  observed  at  VAFB  and  affects  the  heat  energy  balance  at 
ground  level. 

6.  Continue  to  investigate  the  feasibility  of  incorporating  positive  and  negative  source 
buoyancy  effects  in  HOTMAC  and  RAPTAD,  and  test  the  scheme  with  observations. 
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SECTION  m 


WORK  PERFORMED 

A.  TASK  1:  Develop  a  method  to  incorporate  National  Meteorological  Center’s  (NMC) 
Limited-Area  Fine  Mesh  (LFM)  and/or  Nested  Grid  Model  (NGM)  weather 
forecast  into  HOTMAC. 

Accuracy  of  HOTMAC  and  RAPTAD  prediction  should  be  greatly  improved  when  large-scale 
wind  variations  in  time  and  space  are  incorporated  into  HOTMAC  by  using  NMC’s  forecast 
models  LFM  or  NGM. 

Since  the  horizontal  scale  of  the  study  area  (50  kilometers)  is  comparable  to  the  horizontal 
grid  spacing  used  in  the  NMC  LFM  model  (100  kilometers)  or  NGM  (30  kilometers),  large-scale 
wind  distribution  used  to  drive  HOTMAC  will  be  considered  to  be  spatially  uniform  across  the 
horizontal  domain,  though  variable  in  the  vertical  direction. 

Our  recent  experience  (Reference  14)  indicated  that  incorporation  of  large-scale  wind  vari¬ 
ations  (represented  by  upper-air  wind-sounding  data)  into  HOTMAC  improved  simulations  of 
wind  speed  and  wind  profiles  over  complex  terrain  in  the  Sacramento  River  Valley  north  of 
Sacramento,  CA.  We  used  a  method  referred  to  as  Four-Dimensional  Data  Assimilation  (4DDA), 
developed  by  Anthes  (Reference  15)  and  Hoke  and  Anthes  (Reference  16;  Reference  17)  to  in¬ 
corporate  the  large-scale  wind  variations  into  HOTMAC. 

Real-time  forecast  data  were  not  used  in  this  proposal  because  YSA  (Yamada  Science  and 
Ait  Corporation)  could  not  receive  NMC  data  in  real-time.  Thus,  only  archived  data  were  used 
for  development  of  the  method. 

The  following  sections  present  a  brief  description  of  HOTMAC,  simulations  of  wind  distri¬ 
butions  from  HOTMAC,  and  evaluations  of  model  performance  using  Phase  I  WIND  (Wind  In 
Nonuniform  Domain)  data  (Reference  18).  Text  was  duplicated  from  a  manuscnpt  which  was 
submitted  for  publication  to  the  Monthly  Weather  Review  (Reference  14). 
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1 .  Model 


The  basic  equations  for  HOTMAC  are  the  conservation  equations  for  mass,  momentum, 
potential  temperature,  mixing  ratio  of  water  vapor,  and  turbulence  kinetic  energy  (Reference  19). 

The  potential  temperature  equation  was  modified  (Reference  20)  so  that  the  deviation  of 
potential  temperature  from  that  of  the  large-scale  flow  at  an  initial  state  was  solved.  This 
modification  was  necessary  to  maintain  stable  numerical  simulations  and  realistically-predicted 
wind  fields  when  HOTMAC  was  applied  to  simulate  air  flows  over  complex  terrain  with  strong 
wind  shear  and  temperature  inversion  (Reference  19).  The  large-scale  temperature  was  allowed 
to  vary  with  height,  and  assumed  to  be  uniform  in  the  horizontal  directions. 

Also  referred  to  as  a  “second-moment  turbulence-closure  model,”  HOTMAC  is  based  on 
a  set  of  second-moment  turbulence  equations  closed  by  assuming  certain  relationships  between 
unknown  higher-order  turbulence  moments  and  known  lower-order  moments.  HOTMAC  can  be 
used  under  quite  general  conditions  of  flow  and  thermal  stratification:  its  methods  for  turbulence 
parameterization  are  more  advanced  than  those  of  simple  eddy  viscosity  models.  The  present 
model,  which  is  referred  to  as  the  level  2.5  model  (Reference  21),  solves  a  prognostic  equation 
for  turbulence  kinetic  energy  only;  the  remaining  second-moment  turbulence  variables,  such  as 
standard  deviations  of  wind  components,  and  heat  and  momentum  fluxes,  are  solved  from  a  set 
of  algebraic  equations. 

The  present  model  assumes  hydrostatic  equilibrium  and  uses  the  Boussinesq  approximation. 
Therefore,  in  theory,  the  model  applications  are  limited  to  flows  where  the  local  acceleration 
and  advection  terms  in  the  equation  of  vertical  motion  are  much  smaller  than  the  acceleration 
due  to  gravity  (hydrostatic  equilibrium),  and  temperature  variations  in  the  horizontal  directions 
are  not  too  large  (Boussinesq  approximation).  This  assumption  is  probably  justified  given  the 
horizontal  grid  spacing  of  5  kilometers  used  in  this  study.  The  only  way  to  thoroughly  validate 
these  assumptions  would  be  to  repeat  the  simulations  with  a  nonhydrostatic,  non-Boussinesq 
mesoscale  model  and  compare  those  results  with  the  present  results.  Construction  of  a  non- 
Boussinesq  turbulence-closure  model  is  complicated  if  not  impossible. 
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Variations  of  large-scale  wind  distributions  were  incorporated  into  the  equations  of  motion 
through  a  technique  referred  to  as  “nudging”  or  “Newtonian  relaxation  (Reference  15,  Reference 
16).  The  terms  Cn(Ut-U)  and  Cn(VrV)  were  added,  respectively,  to  the  equations  of  motion 
for  the  east- west  and  north-south  components.  Cn  is  the  nudging  coefficient,  and  for  the  present 
study  a  constant  value  of  5  x  10“4  was  employed  throughout  the  simulations.  The  three  levels 
near  the  surface  were  not  subjected  to  nudging.  Ut  and  Vt  are  “target”  wind  components  for  the 
corresponding  wind  components  U  and  V.  Ut  and  Vt  were  computed,  as  in  Equations  (1)  and 
(2),  from  the  large-scale  wind  distributions  and  the  geostrophic  wind,  assuming  a  horizontally 
homogeneous  condition  (Reference  20): 

Ut  =  Uobs  -  r(Vob*  -  Vg)  (1) 

^ n 

Vt  =  V0bs  +  -~{U0bs  ~  Ug )  (2) 

Here,  U0bs  and  Vobs  are  the  observed  wind  components,  and  Ug  and  Vg  the  geostrophic  wind 
components.  The  formulas  for  Ug  and  Vg  can  be  found  in  Yamada  (Reference  22)  and  Yamada 
and  Bunker  (Reference  20).  Ut  and  Vt  are,  in  general,  different  from  the  corresponding  large- 
scale  wind  components.  Observed  winds  may  be  used  as  target  winds  if  the  Coriolis  force  is 
absent  or  if  the  observed  winds  are  identical  to  the  geostrophic  winds.  In  all  other  cases,  if 
only  the  observed  winds  are  used  in  the  nudging,  the  solutions  will  generally  be  diffeient  ftom 
the  observations. 

The  physical  meaning  of  the  target  winds  from  Equations  (1)  and  (2)  is  that  the  solutions 
of  the  equations  of  motion  with  the  target  winds  become  identical  to  the  observed  winds  in  the 
absence  of  frictional  effects.  Thus,  modeled  winds  should  correspond  closely  to  observed  winds 
in  the  layers  above  the  boundary  layer,  where  frictional  effects  are  negligible.  In  the  boundary 
layer  frictional  effects  are  significant,  due  to  atmospheric  turbulence,  and  the  nudging  terms 
play  relatively  minor  roles.  To  ensure  that  the  model  simulates  mesoscale  features,  the  nudging 
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terms  were  turned  off  at  the  first  three  vertical  levels.  In  summary,  the  nudging  terms  force 
the  modeled  winds  to  match  with  observations  in  the  free  atmosphere,  but  they  play  relatively 
minor  roles  in  the  boundary  layer. 

Three-dimensional  fields  of  the  target  wind  were  generated  as  follows:  First,  U0bs  and  V0bS 
were  calculated  from  a  sounding  which  was  interpolated  spatially  to  satisfy  mass  continuity. 
Then,  Uobs  and  Vobs,  together  with  Ug  and  Vg  (computed  by  the  formulas  in  Yamada  (Reference 
22))  were  substituted  into  Equations  (1)  and  (2)  to  obtain  Ut  and  Vt. 

Comparisons  in  the  simulated  wind  fields  were  made  by  nudging  to  the  target  wind 
components  Ut  and  Vt  and  by  nudging  to  the  observed  wind  components  U0bS  and  Vobs.  Nudging 
to  the  target  wind  components  produced  better  agreement  between  simulated  and  observed  upper 
winds  than  nudging  to  the  observed  wind  components.  Therefore,  in  the  following  section,  the 
results  obtained  by  nudging  to  the  target  wind  components  are  shown. 

Surface  boundary  conditions  were  constructed  from  the  empirical  formulas  of  Dyer  and 
Hicks  (Reference  23)  for  nondimensional  wind  and  temperature  profiles.  The  temperatures 
in  the  soil  layers  were  obtained  by  solutions  of  the  heat-conduction  equation.  Appropriate 
boundary  conditions  were  the  heat  balance  at  the  soil  surface  and  specification  of  the  soil 
temperature  at  a  certain  depth.  The  lateral  boundary  values  were  obtained  by  integration  of 
the  corresponding  governing  equations,  except  that  variations  in  the  horizontal  directions  were 
neglected.  Parameterization  of  tall  canopy  effects  on  wind  and  radiation  has  been  studied 
(Reference  24).  These  effects  are  included  in  the  present  model. 

The  governing  equations  were  integrated  by  use  of  the  Alternating  Direction  Implicit 
method  (Reference  25).  A  time  increment  was  chosen  to  be  90  percent  of  the  minimum  value  of 
Axj/Ui,  where  Axj  is  a  grid  spacing  and  Uj  the  velocity  component  in  the  i-th  direction  (Courant- 
Freidrich-Lewy  criterion).  The  integration  time  increment  is  also  limited  by  the  propagation 
speed  of  gravity  waves,  which  is  computed,  based  on  the  potential  temperature  gradients.  To 
increase  the  accuracy  of  finite-difference  approximations,  mean  and  turbulence  variables  are 
defined  at  grids  which  are  staggered  in  both  the  horizontal  and  vertical  directions.  Mean  winds, 
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temperature,  and  water  vapor  vary  most  with  height  near  the  surface.  To  resolve  these  variations 
without  introducing  an  excessive  computational  burden,  nonuniform  grid  spacing  is  used  in  the 
vertical  direction. 

2.  Evaluation  of  Model  Performance 

Visual  comparisons  of  simulation  with  observation  are  useful.  However,  quantitative 
evaluation  of  model  performance  is  also  desirable,  since  it  enables  us  to  objectively  compare 
one  model  to  other  models,  or  one  simulation  case  to  another.  For  the  present  study,  continuous 
observations  of  wind  direction  and  speed,  temperature,  and  humidity  were  available  throughout 
the  24-hour  period  at  21  surface  stations  for  both  the  Phase  I  and  II  measurement  penods,  and 
upper-air  sounding  data  were  available  every  2  hours  at  5  sites. 

Hourly  means,  standard  deviations,  systematic  and  unsystematic  components  of  the  root- 
mean-square  differences  (as  well  as  the  total  root-mean-square  difference),  and  agreement 
measures,  were  calculated  for  surface  wind,  temperature,  and  dew  point  data.  No  upper-air 
statistical  comparison  was  made  due  to  an  insufficient  amount  of  such  data.  Definitions  of  the 
statistical  parameters  are  given  below. 

a.  Mean 


_  Ijizx(k) 

x  =  -^r  (3) 

where  x(k)  is  one  of  the  meteorological  variables  at  the  £-th  station,  and  N  is  the  total  numbei  of 
stations.  Means  for  both  the  simulation  and  the  observation  were  calculated;  for  good  agreement, 
both  should  be  similar. 

b.  Standard  deviation 


a  = 


[xjk)  -  xf  }  ' 
“  N  1 


(4) 
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The  standard  deviation  for  both  the  simulation  and  the  observation  should  also  be  similar  for 
good  agreement. 

c.  Root-mean-square  differences 


Use  of  systematic  and  unsystematic  components  of  the  root-mean-square  differences  as 
well  as  the  total  root-mean-square  difference  (RMSDs,  RMSDU,  and  RMSD,  respectively)  was 
recommended  by  Willmott  (Reference  26)  for  quantitative  evaluation  of  model  performance. 
Steyn  and  McKendry  (Reference  27)  used  these  parameters  to  evaluate  the  performance  of  the 
Colorado  State  University  mesoscale  model  by  Mahrer  and  Pie  Ike  (Reference  28;  Reference  29) 
for  a  sea-breeze  case  over  complex  terrain.  These  parameters  were  used  also  by  Ulrickson  and 
Mass  (Reference  30)  for  the  evaluation  of  the  above-mentioned  model  when  used  for  simulation 
of  mesoscale  circulations  over  the  Los  Angeles  Basin. 

These  parameters  are  defined  as: 

1 
2 

(5) 


RMSDs  - 


is -  x0(k)f 


RMSDU 


N‘ 


Zk(x*m(k)  -  xm(k)Y 


(6) 


RMSD  = 


—T,k(xm(k)  x0(k)) 


(7). 


where  N  is  the  number  of  stations  (evaluation  points),  xm  and  x0  are  simulated  and  observed 
values  respectively,  and  x*m  =  a  +  bx0  where  a  and  b  are  the  parameters  associated  with  an 
ordinary  least  squares  linear  regression  between  x0  and  x,„  (Reference  27). 

RMSDs  is  an  estimate  not  only  of  the  offset  bias  of  a  model,  but  also  of  any  linear 
variation  in  the  bias  of  that  model.  RMSDU  is  a  measure  of  the  nonlinear  discrepancy  between 
simulation  and  observations  (Reference  30).  The  definitions  of  systematic  and  unsystematic 
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RMSD  apply  to  the  spatial  coherence  of  each  hourly  set  of  simulations  and  observations,  and 
do  not  address  the  errors  at  individual  sites. 

d.  Agreement  measure 


Sjfc( |:rm( fc)  %m\  T  \x0(k)  ^o|) 

This  dimensionless  index  has  a  theoretical  range  of  1.0  (for  perfect  agreement)  to  0.0  (for  no 
agreement).  This  parameter  was  also  used  by  Steyn  and  McKendry  (Reference  27),  and  Ulrickson 
and  Mass  (Reference  30)  for  their  studies  of  model  evaluation. 

3.  Simulation  of  Phase  I  Data 

HOTMAC  was  used  for  simulation  of  the  project  WIND  Phase  I  (summer)  data  (Reference 
18).  Phase  I  of  the  field  study  was  conducted  in  June  and  July  1985  during  a  weak-to-moderate 
marine  incursion  regime.  The  main  meteorological  events  during  the  Phase  I  observation  period 
were  the  diurnal  variation  of  the  surface  wind  fields,  upslope  winds  during  the  day  and  downslope 
winds  during  the  night.  Simulation  of  this  diurnal  change  —  its  timing,  strength,  the  structure 
of  the  winds,  and  the  accompanying  changes  in  temperature  near  the  surface  —  was  the  mam 
focus  of  this  study.1 

Simulation  results  were  compared  temporally  and  spatially  with  both  surface  and  upper-air 
observations.  Statistical  comparisons  between  simulation  and  observation  were  also  made.  The 
results  indicated  good  simulation  capabilities  of  the  model. 

a.  Model  Domain 

The  study  area  was  centered  in  the  Sacramento  River  Valley  north  of  Sacramento,  CA, 
extending  eastward  up  the  slopes  of  the  Sierra  Nevada  and  westward  to  the  slopes  of  the  Coast 
Range.  Latitude  and  longitude  of  the  southwestern  comer  of  the  domain  for  the  simulation 

1  Pearce.  R.  (Personal  communication),  1992. 


18 


(Figure  4)  were  39.63  N  and  122.99  W.  Meteorological  variables  were  calculated  at  40  X  40 
grid  points  with  a  horizontal  grid  spacing  of  5  kilometers.  In  the  figure,  characters  represent 
the  locations  of  surface  stations.  Upper-air  soundings  were  taken  near  BC3(01),  C7(02),  S8(03), 
B2(04),  and  S3(05).  Here,  numbers  in  parentheses  represent  upper-air  stations. 

b.  Initial  Values 

A  meteorological  data  set  for  the  24-hour  period  starting  at  0900  LT  of  Julian  day 
178  (27  June)  was  used  in  this  study,  this  day  being  characterized  as  a  typical  summer  day. 
The  most  common  flow  regime  in  the  Sacramento  River  Valley  during  the  summer  is  southerly 
(upslope)  flow  during  the  day  and  northerly  (downslope)  flow  at  night.  The  southerly  flow  is 
consistent  with  a  weak  marine  incursion,  which  results  when  wind  that  blows  onshore  from  the 
San  Francisco  Bay  area  (more  than  200  kilometers  away)  moves  inland  and  becomes  diverted 
north  and  south  by  the  Sierra  Nevada  along  the  Sacramento  River  Valley. 

An  initial  wind  profile  at  Station  3  was  first  constructed  by  assuming  a  logarithmic 
variation  with  u*  =  0.2  m/s  and  z0  =  0.1m  from  the  ground  up  to  the  level  where  the  wind  speed 
reaches  an  ambient  value  of  5  m/s.  Initial  wind  profiles  at  other  grid  locations  were  obtained  by 
multiplying  the  station  3  winds  by  (z  -  zg)/(H  -  zg)  to  satisfy  mass  continuity  approximately,  and 
wind  directions  were  initially  from  the  SW  (210  degrees)  throughout  the  atmosphere.  The  vertical 
profile  of  potential  temperature  was  initialized  based  on  an  upper-air  sounding  taken  at  station 
3  at  0900  LT  of  day  178,  and  initial  potential  temperatures  were  assumed  to  be  uniform  in  the 
horizontal  directions.  Initial  values  of  water  vapor  were  also  based  on  observation,  whereas  the 
turbulence  kinetic  energy  and  length  scale  were  initialized  using  the  initial  wind  and  temperature 
profiles,  and  the  relationships  determined  by  the  Level  2  turbulence-closure  model  (Reference 
31).  The  Level  2  model  assumes  a  balance  between  the  production  and  dissipation  terms  in  the 
turbulence  kinetic  energy  equation. 

The  simulation  commenced  at  0900  LT  of  day  178  and  lasted  until  0900  of  day  179. 
However,  model  calculations  were  begun  4  hours  before  the  start  of  the  actual  simulation  to 
produce  an  initial  adjustment.  At  1  hour  before  the  start  of  simulation,  the  observed  wind  data  at 
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Figure  4.  Model  Domain  (200  X  200  Kilometers)  for  WIND  Phase  I  Simulation.  Terrain  is 
Contoured  by  Solid  Lines  in  Increments  of  400  Meters.  The  Lowest  Contour  is  400 
Meters  above  Mean  Sea  Level.  The  Locations  of  Surface  Stations  are  Indicated  by 
Characters.  Numbers  in  Parentheses  represent  Upper-Air  Stations. 
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0900  were  read  in  and  model  winds  began  nudging  toward  target  winds,  as  given  in  Equations 
(1)  and  (2). 

c.  Surface  Meteorological  Parameters 

(1)  Horizontal  Wind  Vector  Fields.  Model  calculations  were  made  over  the  area  of  200 
X  200  kilometers  as  shown  in  Figure  4.  However,  to  make  visual  comparison  easy,  simulated 
horizontal  wind  vectors  were  plotted  only  for  the  area  of  100  X  100  kilometers  that  contains 
all  of  the  surface  and  upper-air  stations.  Figures  5  and  6  show  surface  (10-meter  level)  wind 
vector  fields  at  1500  LT  of  day  178,  and  0500  LT  of  day  179.  In  these  figures,  (A)  represents 
the  simulated  wind  field  with  nudging  to  the  target  winds,  (B)  the  observed  winds,  (C)  the 
target  wind  field,  and  (D)  the  simulated  wind  field  without  nudging.  An  arrow  with  a  unit  grid 
distance  represents  a  wind  vector  of  10  m/sec.  The  four  fields  allow  us  to  examine  the  influences 
of  nudging.  It  can  be  seen  in  both  figures  that  nudging  to  the  target  wind  fields  has  significantly 
influenced  the  simulated  wind  fields.  In  Figure  5,  nudging  to  the  target  wind  field  has  produced 
southerly  flow  throughout  the  model  domain,  whereas  the  simulated  wind  field  without  nudging 
was  southwesterly.  In  Figure  6,  the  effects  of  nudging  to  the  target  wind  field  were  to  shift  the 
wind  direction  from  easterly  along  the  Sierra  Nevada  to  northeasterly,  and  to  generate  northerly 
flow  in  the  Sacramento  River  Valley.  It  was  also  noted,  from  the  comparisons  between  (A)  and 
(D),  that  the  simulations  with  nudging  to  the  target  wind  fields  produced  smoother  wind  fields 
throughout  the  domain  than  those  without  nudging. 

Upslope  winds  during  the  day  (Figure  5(A))  and  downslope  winds  dunng  the  night 
(Figure  6(A))  can  be  detected  along  the  foothills  of  the  Sierra  Nevada.  The  observations  show 
that  during  the  day,  wind  vector  distributions  were  relatively  uniform,  but  during  the  night, 
wind  vectors  varied  considerably  in  space,  because  the  surface  layer  was  decoupled  from  the 
atmosphere  above,  and  locally  generated  forces  dominated  the  wind  distributions.  These  local 
forces  appeared  to  be  due  to  subgrid  scale  phenomena.  Neither  simulations  with  nor  without 
nudging  produced  the  large  variations  of  wind  distributions  detected  in  the  observations. 

(2)  Time  Series  Comparisons  of  Station  Data.  The  evolution  of  the  wind  direction  and 
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Figure  5.  Surface  (10-Meter  Level)  Horizontal  Wind  Vector  Distributions  for  Phase  I  at  1500 
LT  of  Day  178.  (A)  the  Simulation  With  Nudging,  (B)  Observations,  (C)  the  Target 
Wind  Field,  and  (D)  the  Simulation  Without  Nudging.  An  Arrow  With  a  Unit  Gnd 
Distance  Represents  10  m/sec. 
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Figure  6.  Same  as  Figure  5  except  at  0500  LT  of  Day  179. 
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speed  at  the  10-meter  level,  and  of  the  temperature  and  dew  point  at  the  2-meter  level  are  given 
in  Figures  7(A),  (B),  (C),  and  (D)  for  four  representative  stations,  S2,  S10,  S14,  and  C3.  These 
figures  show  the  results  of  simulation  with  nudging  to  the  target  wind.  In  these  figures,  hourly 
values  of  observations  are  plotted  with  circles,  and  simulations  by  solid  lines. 

Station  S2  was  located  in  the  valley.  The  wind  direction  shift,  from  south  early  during 
the  daytime  to  the  north  early  during  the  nighttime,  was  simulated  well.  The  diurnal  variation 
of  the  modeled  wind  speed  is  also  in  good  agreement  with  the  observations. 

Station  S10,  which  was  located  near  the  foothills  of  the  Sienra  Nevada,  showed 
disagreement  between  simulated  and  observed  wind  directions  after  sunset.  Whereas  the  observed 
wind  direction  was  predominantly  northwesterly,  the  simulated  direction  was  northeasterly. 
Observed  wind  speeds  during  the  nocturnal  period  were  approximately  1  m/s,  and  at  such  low 
speeds  determination  of  wind  direction  is  difficult  and  may  contain  large  errors. 

At  Station  SI 4,  which  was  located  on  the  western  slope  of  the  mountain  range,  the 
observed  wind  direction  changed  gradually  from  southeast  to  east,  but  the  simulated  direction 
showed  a  sudden  shift  of  wind  direction  from  south  to  east  after  sunset.  Wind  speeds  did  not 
agree  well  during  the  day. 

At  Station  C3,  the  observed  wind  direction  changed  from  southwest  early  in  the  day  to 
northwest  early  at  night,  whereas  the  simulated  wind  direction  stayed  southerly  throughout  the 
24-hour  period.  The  wind  speeds  also  did  not  agree  well.  At  this  station,  the  observed  wind 
speeds  stayed  very  low  during  the  night;  the  model  failed  to  simulate  these  low  wind  speeds. 

The  simulated  temperatures  agreed  well  with  observations  during  the  day,  but  during 
the  night  the  differences  between  simulated  and  observed  temperatures  became  significant. 
Particularly  at  C3,  the  model  failed  to  simulate  the  nocturnal  temperature  drop.  Station  C3 
was  located  in  a  meadow  which  was  surrounded  by  high  mountains.  This  is  clearly  an  example 
of  subgrid  scale  phenomena,  which  the  model  could  not  resolve  with  a  5-kilometer  gnd  spacing. 
In  general,  the  diurnal  ranges  of  simulated  temperature  were  smaller  than  those  observed. 

Throughout  the  24-hour  period,  the  dew  points  remained  low,  as  can  be  seen  from 
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Figure  7.  Time  Evolutions  of  Surface  Meteorological  Variables  During  the  24-Hour  Period  of 
the  Project  WIND  Phase  I,  at  (A)  Station  S2,  (B)  Station  S10,  (C)  Station  S14,  and 
(D)  Station  C3.  (-  Simulation,  and  o  Observation)  (Continued) 
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Figure  7.  Time  Evolutions  of  Surface  Meteorological  Variables  During  the  24-Hour  Period  of 
the  Project  WIND  Phase  I,  at  (A)  Station  S2,  (B)  Station  S10,  (C)  Station  S14,  and 
(D)  Station  C3.  (-  Simulation,  and  o  Observation)  (Concluded) 
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Figures  7(A),  (B),  and  (C).  There  was  no  dew-point  observation  available  at  C3. 

(3)  Statistical  Parameters.  Using  the  data  from  2 1  surface  observation  sites,  mean  wind 
directions  and  speeds,  in  addition  to  standard  deviation  of  speed,  were  calculated  hourly  for  the 
simulation  (solid  line)  and  observation  (circles)  as  shown  in  Figure  8(A).  Also  shown  are  hourly 
RMSD,  RMSDs  and  RMSDU,  and  agreement  measures  of  wind  speeds.  Mean  wind  direction 
was  calculated  from  the  means  of  the  horizontal  wind  components. 

The  top  two  portions  of  Figure  8(A)  show  that  the  overall  patterns  of  the  observed  and 
simulated  surface  wind  fields  agreed  well,  except  for  several  hours  after  sunset. 

Mean  directions  of  observed  wind  were  from  the  southwest  during  the  day,  shifted 
gradually  to  the  north  after  sunset,  and  stayed  northeast  during  the  night.  The  simulated  mean 
wind  directions  agreed  well  with  the  observed  directions  during  the  day,  but  the  shifting  to  the 
north  was  delayed  for  several  hours  in  the  simulation,  and  after  midnight  stayed  north  rather  than 
northeast.  Simulated  and  observed  mean  wind  speeds  agreed  well  throughout  the  24-ho.ur  period. 

The  standard  deviations  of  the  observed  and  simulated  wind  speeds  were  comparable 
throughout  the  simulation  period.  The  evolution  of  RMSD  and  its  two  components  shows  no 
diurnal  trend,  and  both  components  of  RMSD  are  of  roughly  equal  magnitude.  The  unsystematic 
component,  RMSDU,  represents  the  irreducible  deviation  between  observed  and  simulated  results, 
while  the  systematic  component  represents  trend- like  differences  between  observed  and  simulated 
fields  (Reference  27). 

The  greatest  agreement  measure  (0.65)  was  obtained  at  1100  LT  of  Day  178,  but  after 
1200  LT  it  remained  near  0.5,  without  any  systematic  trend. 

Similar  statistical  calculations  were  conducted  for  the  simulated  wind  fields  without 
nudging  and  for  the  persistence  wind  field,  and  the  results  are  shown  in  Figure  8(B)  and  (C). 
Here,  the  persistence  wind  field  was  the  observed  10-meter  level  wind  field  at  0900  LT  of  Day 
178,  when  the  simulation  was  initialized.  From  the  comparisons  between  Figures  8(A),  (B),  and 
(C),  the  following  can  be  inferred: 
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Figure  8  Time  Evolutions  of  Statistical  Parameters  of  Surface  (10-Meter  Level)  Wind  Speed 
and  Direction  for  24-Hour  Period  of  Phase  I  for  the  Simulation  with  Nudging  to  the 
Target  Winds.  From  the  Top,  Mean  Wind  Direction  (  -  Simulation,  and  o  Observation), 
Mean  Wind  Speed  (  -  Simulation,  and  o  Observation),  Standard  Deviation  of  Wind 
Speed  (  -  Simulation,  and  o  Observation),  Root-Mean-Square  Differences  (*  RMSD, 
x  RMSDu,  and  o  RMSDs),  and  Agreement  Measure  of  Wind  Speed.  (A)  Winds  were 
Nudged  to  the  Target  Winds,  (B)  No  Nudging  was  applied,  and  (C)  Persistence  Winds 
were  used.  (Continued) 
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Figure  8.  Time  Evolutions  of  Statistical  Parameters  of  Surface  (10-Meter  Level)  Wind  Speed 
and  Direction  for  24-Hour  Period  of  Phase  I  for  the  Simulation  with  Nudging  to  the 
Target  Winds.  From  the  Top,  Mean  Wind  Direction  (-  Simulation,  and  o  Observation), 
Mean  Wind  Speed  (-  Simulation,  and  o  Observation),  Standard  Deviation  of  Wind 
Speed  (-  Simulation,  and  o  Observation),  Root-Mean-Square  Differences  (*  RMSD, 
x  RMSDU,  and  o  RMSDs),  and  Agreement  Measure  of  Wind  Speed.  (A)  Winds  were 
Nudged  to  the  Target  Winds,  (B)  No  Nudging  was  applied,  and  (C)  Persistence  Winds 
were  used.  (Continued) 
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Figure  8.  Time  Evolutions  of  Statistical  Parameters  of  Surface  (10-Meter  Level)  Wind  Speed 
and  Direction  for  24-Hour  Period  of  Phase  I  for  the  Simulation  with  Nudging  to  the 
Target  Winds.  From  the  Top,  Mean  Wind  Direction  (-  Simulation,  and  o  Observation), 
Mean  Wind  Speed  (-  Simulation,  and  o  Observation),  Standard  Deviation  of  Wind 
Speed  (-  Simulation,  and  o  Observation),  Root-Mean-Square  Differences  (*  RMSD, 
x  RMS Du,  and  o  RMSDs),  and  Agreement  Measure  of  Wind  Speed.  (A)  Winds  were 
Nudged  to  the  Target  Winds,  (B)  No  Nudging  was  applied,  and  (C)  Persistence  Winds 
were  used.  (Concluded) 
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(1)  The  simulation  with  nudging  to  the  target  winds  resulted  in  better  agreements  in  mean 
wind  directions  and  speeds  than  the  simulation  without  nudging  or  the  persistence. 

(2)  The  standard  deviations  and  the  RMSDs  of  the  three  different  wind  fields  were  of  similar 
magnitudes. 

(3)  The  agreement  measures  of  wind  speed  calculated  for  the  persistence  wind  were  the 
greatest  among  the  three  different  wind  fields. 

Thus,  an  evaluation  of  model  performance  should  not  be  based  on  a  single  statistical 
parameter.  The  use  of  an  increased  number  of  observation  sites,  distributed  evenly  over  the 
model  domain,  would  yield  a  better  statistical  evaluation. 

Mean,  standard  deviation,  RMSD  and  agreement  measure  of  temperature  at  the  2-meter 
level  were  also  calculated  hourly  as  shown  in  Figure  9.  Good  agreement  existed  between 
the  simulations  and  observations  during  the  day,  but  during  the  night  there  were  noticeable 
differences  between  the  observed  and  simulated  fields.  The  diurnal  variation  of  the  simulated 
surface  temperature  was  much  smaller  than  that  of  the  observed,  resulting  in  poorer  agreement 
during  the  night  than  during  the  day.  Smaller  variations  in  the  simulated  temperature  field 
than  those  of  the  observed  temperature  field  resulted  in  large  values  of  RMSDs  during  the 
nocturnal  period.  Surface  temperature  is  dependent  on  many  physical  parameters,  including 
surface  vegetation  cover  and  the  thermal  conductivity  of  the  soil.  However,  in  the  present  study, 
a  bare  surface  condition  and  a  constant  thermal  conductivity  were  assumed  throughout  the  model 
domain.  The  incorporation  of  locally-dependent  surface  characteristics  might  have  improved  the 
prediction  of  surface  temperature. 

Temperature  fields  of  the  simulation  without  nudging  yielded  values  of  the  statistical 
parameters  very  similar  to  those  with  nudging  (the  results  are  not  shown). 

d.  Upper-Air  Variables 

To  compare  upper-air  variables  at  certain  levels,  the  simulated  and  observed  vertical 
profiles  were  expressed  by  bi-cubic  spline  curves.  Figure  10  shows  the  evolution  of 
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Figure  9.  Time  Evolutions  of  Statistical  Parameters  of  Surface  (2-Meter  Level)  Temperature. 
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Simulation,  and  o  Observation),  Root-Mean-Square  Differences  (*RMSD,  x  RMSDU, 
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Figure  10.  Time  Evolution  of  Wind  Direction,  Wind  Speed,  Temperature  and  Dew  Points  at 
Different  Heights  (agl)  at  Upper-Air  Station  1  during  the  24-Hour  period  of  Phase 
I  (-  Simulation,  and  *  Observation). 
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upper-air  variables  (wind  direction  and  wind  speed,  temperature,  and  dew  point)  at  six  different 
height  levels  (4000,  2000,  1000,  200,  50,  and  10  meters  agl)  for  Station  1.  Note  that  the  upper- 
air  data  taken  at  Station  3  were  assimilated  into  the  model,  and  that  Station  1  was  located  about 
30  kilometers  east  of  Station  3  on  the  western  slope  of  the  Sierra  Nevada. 

The  agreement  between  model  and  observation  was  excellent  throughout  the  simulation 
at  all  levels,  as  seen  in  these  figures,  and  the  effects  of  diurnal  heating  and  cooling  were  clearly 
recognized  at  the  lower  levels.  Although  temperatures  and  dew  points  in  the  layers  above 
200  meters  stayed  almost  constant  throughout  the  24-hour  period,  wind  directions  and  speeds 
showed  some  variation. 

Vertical  distributions  of  wind  direction  and  wind  speed,  temperature,  and  dew  point  at 
1 100  and  1700  LT  of  Day  178,  and  0100  and  0700  LT  of  Day  179  at  Station  4  are  shown  in  Figure 
11.  In  Figure  11,  simulations  are  plotted  by  solid  lines  and  observations  by  circles.  Excellent 
agreement  in  vertical  distributions  of  wind  speed  and  direction,  temperature  and  dew  point  were 
obtained  throughout  the  simulation,  and  good  agreements  were  also  obtained  at  other  stations. 
As  noted  previously,  no  nudging  was  performed  to  the  temperature  and  mixing  ratio  equations. 

e.  Summary 

A  three-dimensional  mesoscale  model,  HOTMAC,  was  used  to  simulate  the  project 
WIND  Phase  I  data.  The  nudging  method  in  which  horizontal  wind  components  were  nudged 
to  the  target  wind  components  defined  by  Equations  (1)  and  (2)  was  used.  The  advantage  of 
nudging  to  the  target  winds  over  nudging  to  the  observed  winds  was  illustrated. 

Phase  I  was  characterized  by  typical  summer  anticyclonic  conditions  over  the  Sacra¬ 
mento  River  Valley  in  northern  California.  Southerly  flow  (upslope  wind)  during  the  day  and 
northerly  flow  (downslope  wind)  at  night  were  the  characteristic  flow  patterns. 

Comparisons  between  observations  and  simulations  of  the  wind  duections  and  speed, 
temperature,  and  dew  point  showed  fairly  good  agreement,  both  spatially  and  temporally, 
throughout  entire  24-hour  periods.  Thus,  a  single  station  sounding  could  be  used  to  represent 
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Figure  11.  Vertical  Distributions  of  Wind  Direction  and  Speed,  Temperature  and  Dew  Point  at 
Station  4  for  Phase  I.  (A)  1100  LT,  Day  178,  (B)  1700  LT,  Day  178,  (C)  0100  LT, 
Day  179,  (D)  0700  LT,  Day  179  (-  Simulation,  and  o  Obseiwation). 


large-scale  meteorological  conditions  over  the  model  domain  fairly  well.  However,  if  the  area 
were  under  strong  synoptic  influences,  the  present  method  should  not  be  applied,  and  synoptic- 
scale  weather  variations  should  be  incorporated  into  a  mesoscale  model  by  specifying  time- 
dependent  lateral  boundary  conditions  or  by  nesting  a  mesoscale  model  into  a  large-scale  model. 

Statistical  parameters,  including  the  mean  and  standard  deviation  of  both  the  observed 
and  simulated  variables,  the  systematic  and  unsystematic  components  of  the  root-mean-square 
differences  as  well  as  the  total  root-mean- square  difference  (RMSDs,  RMSDu,  and  RMSD, 
respectively),  and  the  agreement  measures  were  calculated  hourly  using  the  surface  station  data. 
These  statistical  parameters  were  calculated  for  three  different  wind  fields  simulated  by  nudging 
to  the  target  winds,  with  no  nudging,  and  with  persistence.  The  wind  fields  obtained  by  nudging  to 
target  winds  yielded  better  agreement  in  wind  direction  and  speed  when  compared  to  the  two  other 
wind  fields.  However,  the  magnitudes  of  the  standard  deviation  and  RMSD  were  comparable 
among  the  three  wind  fields.  Agreement  measures  of  wind  speed  using  the  persistence  winds 
were  greater  than  those  obtained  from  the  two  other  wind  fields.  Hence,  model  performance 
should  be  evaluated  carefully,  by  several  different  means. 

B.  TASK  2:  Develop  a  method  to  restart  HOTMAC  prediction  and  maintain  on  disk  the  wind 
and  turbulence  forecast  data  for  the  next  24  hours. 

A  complex  scheme  was  developed  to  operate  HOTMAC  and  RAPTAD  interactively.  HOT¬ 
MAC  runs  continuously  and  always  maintains  24-hour  forecast  data.  The  interactive  scheme 
dictates  that  as  soon  as  RAPTAD  starts  computation,  HOTMAC  halts  its  forecasting.  RAPTAD 
displays  an  hourly  surface-concentration  distribution  plot  on  the  screen  as  soon  as  one  hourly 
computation  is  completed. 

HOTMAC  produces  24  hourly  files,  hotOl  through  hot24.  Once  a  24-hour  forecast  is 
completed  and  the  forecast  data  are  written  to  a  disk,  HOTMAC  “sleeps”  until  the  next 
computation  time.  HOTMAC  “wakes  up”  at  the  clock  hour  and  begins  to  advance  the  forecast 
one  hour.  At  the  end  of  the  one  hour  forecast,  files  are  renamed:  hot02  becomes  hotOl,  hot03 
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becomes  hot02,  and  so  on.  The  newest  forecast  is  written  as  hot24,  so  that  forecast  data  for 
the  next  24  hours  are  always  available. 

RAPTAD  reads  HOTMAC  output  files  hotxx,  where  wind  speed,  wind  direction,  temperature, 
and  turbulence  forecast  data  are  stored.  Graphics  program  hotplt  produces  graphics  for  vertical 
profiles  of  wind,  temperature,  and  turbulence. 

If  the  continuous  operation  of  HOTMAC  is  interrupted  (for  maintenance,  or  by  a  power  failure 
or  mistakes),  HOTMAC  may  be  restarted  by  using  one  of  the  hotxx  files.  The  restart  time  may 
be  specified  in  input  file  hotinp2.  The  restart  parameter  may  be  set  to  zero  if  reinitialization 
is  desired. 

In  addition  to  using  the  default  mode’s  24-hour  forecast,  the  user  can  extend  the  range  of 
the  forecast  period,  up  to  48  hours,  for  planning  operations. 

C.  TASK  3:  Develop  a  method  to  combine  upper-air  wind-sounding  and  tower  data  taken  at 
VAFB  with  HOTMAC  forecast  winds  so  that  RAPTAD  can  use  the  composite 
(forecast  and  observed)  winds. 

Winds  measured  by  upper-air  soundings  and  at  towers  were  combined  with  the  winds 
forecasted  by  HOTMAC.  A  simple  i/r2  weighting  method  was  used  to  compute  the  wind  at 
a  puff  center.  First,  we  computed  wind  at  a  puff  center  based  on  the  winds  measured  at  towers. 
For  example.  Figures  12  through  16  show  wind  vectors  at  tower  locations  (with  numbers) 
and  those  extrapolated  for  puff  centers  (without  numbers)  at  different  hours.  Tower  winds 
show  considerable  variation  in  space,  but  winds  at  puff  centers  are  relatively  uniform,  because 
variations  of  tower  winds  are  minimized  by  averaging.  These  wind  distributions  are  similar 
to  those  generated  by  diagnostic  wind  models  where  winds  are  adjusted  to  satisfy  the  mass 
continuity  equation. 

The  winds  at  puff  centers  (obtained  from  tower  winds)  were  then  combined  with  winds 
computed  by  HOTMAC.  Figures  17-20  show  examples  of  composed  wind  vectors  at  puff  centers. 
Tower  winds  are  again  indicated  by  wind  vectors  with  numbers.  Tower  winds  are  weighted  based 


37 


y  [km] 


WIND  VECTORS 
VAFB  1 1 393  1110  12ft 


x  [km] 


Figure  14.  Similar  to  Figure  12  but  for  January  13,  1993. 
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Figure  15.  Similar  to  Figure  12  but  for  January  14,  1993. 
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Figure  17.  Wind  Vectors  at  Tower  Locations  (With  Numbers)  and  Those  at  Puff  Centers  (Without 
Numbers)  at  0300  Z,  January  12,  1993.  Puff  Velocities  Were  Computed  from  Tower 
Data  and  HOTMAC  Winds  by  1/r2  Weighting  Method. 
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Figure  19.  Similar  to  Figure  17  but  for  0900  Z,  January  12,  1993. 
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on  the  distance  between  a  puff  center  and  the  nearest  tower  site.  If  a  puff  is  located  at  a  tower  site, 
then  the  tower  wind  is  used  for  puff  velocity.  The  weight  function  decreases  exponentially  as 
the  distance  between  a  puff  and  a  tower  increases.  Such  an  assumption  is  particularly  reasonable 
for  airflows  over  complex  terrain,  where  wind  distributions  are  largely  localized. 

Tower  data  are  optional.  An  input  file  rpinp  sets  a  flag  for  simulations  with  and  without 
tower  data. 

D.  TASK  4:  Add  to  the  model  the  physics  necessary  to  forecast  the  formation  and  dissipation 
of  fog. 

Fog  and/or  low  stratus  clouds  are  frequently  observed  at  VAFB,  and  affect  the  heat  energy 
balance  at  the  ground.  The  interaction  of  fog,  solar  heating,  long-wave  radiation  heating/cooling, 
and  turbulence  mixing  is  complex  and  is  not  yet  well  understood.  Fog  physics  was  incorporated 
into  a  one-dimensional  version  of  HOTMAC  and  the  results  were  tested  with  the  data  taken  at 
Cabaw  meteorological  tower  (Reference  32).  The  results  were  satisfactory,  and  the  fog  physics 
was  incorporated  into  a  three-dimensional  version  of  HOTMAC. 

1.  Ensemble  Cloud  Modeling 

The  terms  for  the  rate  of  condensation  are  purposely  eliminated  by  introducing  the  liquid 
water  potential  temperature  (0/)  and  mixing  ratio  of  total  wateT  (Qi).  In  order  to  recover 
the  potential  (or  absolute)  temperature  and  the  mixing  ratios  of  water  vapor  and  liquid  water, 
Gaussian  cloud  relations,  proposed  by  Sommeria  and  Deardorff  (Reference  33)  and  Mellor 
(Reference  34),  are  used.  The  present  method  has  been  applied  to  simulations  of  the  BOMEX 
data  (Reference  35),  and  GATE  data  (Reference  36). 

Following  Mellor  (Reference  35)  we  assume  the  probability  density  function  G  for  0;  and 
Qi  to  be  Gaussian  so  that 
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where  4  =  d2  ,<j2w  =  q2w,  r  =  9iqw/aeiaqWt  and  ~  represents  an  instantaneous  value.  The 
results  are  presented  here  without  details  of  derivation.  Readers  are  referred  to  Mellor  (Reference 
34)  for  details. 

First  we  define 
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where  &/  is  the  mean  saturation  mixing  ratio  of  water  vapor  at  7);  Ts  is  the  liquid  water 
temperature  defined  by  Equation  (16);  es(Ti)  is  the  saturation  water  vapor  pressure  at  7),  and  Rw 
is  the  gas  constant  for  water  vapor.  Mellor  (Reference  34)  obtained  the  following  expiessions.  A 
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function  R  which  indicates  a  fraction  of  cloud  coverage  for  a  given  volume  of  air  is  found,  i.e., 
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Finally,  qf,  the  variance  of  cloud  water  mixing  ratio,  is  expressed  as 


gc 

4(7,  2 


=  R 


l  +  (  Q  i  — 


Qc 

2a, 


+  (V# 


2(7 3 )  y/2jr 


exp 


Ql 1 

2 


(17) 


(18) 


(19) 


(20) 


(21a,b) 


(22) 


(23) 


(24) 


Figure  21  shows  R,  Qc/'2as,  and  q2 /4a2  as  functions  of  Qi  obtained  according  to  Equations 
(17),  (18),  and  (24),  respectively. 
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Figure  21.  Fraction  of  Cloud  Coverage  R,  Mixing  Ration  of  Liquid  Water 
Qc ,  and  Variance  of  Liquid  Water,  as  a  Function  of  Qi. 


As  seen  from  Figure  21,  a  fraction  of  cloud  coverage  R  vanes  gradually  with  Q\  and  takes 
non-zero  values  even  if  Qi  is  negative.  In  other  words,  clouds  can  exist  even  if  the  mixing 
ratio  of  water  vapor  averaged  over  a  gnd  volume  is  not  saturated.  This  is  realistic  because 
the  gnd  spacing  normally  used  in  mesoscale  models  are  larger  than  the  size  of  small  clouds. 
Therefore,  the  present  cloud  model  can  use  a  relatively  large  grid  spacing,  which  could  save 
substantial  computational  time.  A  statistical  cloud  model  such  as  the  present  one  avoids  the 
ambiguous  condensation  cnteria  often  used  by  coarse-gnd  models,  where  saturation  values  are 
lowered  arbitrarily  to  compensate  for  the  amount  of  cloud  not  resolved  by  the  gnds. 

The  turbulence  second  moments  in  Equation  (20)  and  Equations  (21  a,b)  are  obtained  by 
solving  second-moment  turbulence-closure  equations.  Those  expressions  are  given  in  Yamada 
and  Mellor  (Reference  35). 

2.  Short-Wave  and  Long- Wave  Radiation  Parametenzations 

Solar  and  long-wave  radiation  play  important  roles  in  the  formation  and  dissipation  of 
clouds  and  fog.  For  example,  radiation  fog  forms  when  moist  air  near  the  surface  condenses 
as  the  ground  temperature  decreases  due  to  long-wave  radiation  cooling.  As  the  sun  rises,  fog 
reflects,  absorbs,  and  transfers  short-wave  solar  radiation.  The  solar  energy  absorbed  m  a  fog 
layer  heats  the  air  and  converts  liquid  water  to  water  vapor.  The  solar  energy  transmitted  through 
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fog  heats  the  ground  and  increases  ground  temperature.  The  warmed  ground,  in  turn,  heats  the 
air  above  it. 

This  interaction  between  radiative  transfer  and  cloud  development  is  thus  understood  qual¬ 
itatively.  Quantitative  description  is  complex,  however,  and  involves  considerable  computations 
(Reference  37;  Reference  38).  This  is  why  many  numerical  models  either  neglect  the  radiative 
effects  on  clouds  or  adopt  a  very  simple,  highly  parameterized  method. 

Hanson  and  Derr  (Reference  39)  proposed  a  parameterized  solar  radiation  scheme  whose 
parameters  were  obtained  by  curve-fitting  to  numencal  radiative-transfer  results  using  the 
ATRAD  narrow-band  model  (Reference  40).  This  parameterized  method  is  simple,  yet  re¬ 
produced  solar  flux  profiles  within  a  single  cloud  layer  which  were  in  good  agreement  with  the 
numerical  results  obtained  from  ATRAD. 

Motivated  by  the  success  of  the  solar  radiation  scheme,  Hanson  and  Derr  proposed  an 
infrared  (IR)  radiation  scheme  expressed  by  exponential  functions  whose  decay  parameters  were 
determined  by  the  emissivity  methods  (Reference  41).  This  method  reproduced  IR  flux  profiles 
within  layered  clouds  which  were  in  good  agreement  with  the  numerical  results  and  observations 
for  thick  clouds  (-800  meters).  However,  the  parameterization  overestimated  the  flux  decrease 
at  the  cloud  top  and  underestimated  the  cloud  base  warming  compared  to  the  numerical  model 
for  thin  clouds  (-300  meters). 

We  proposed  to  adopt  the  Hanson  and  Derr  parameterization  scheme  because  of  its 
simplicity  and  because  it  can  produce  flux  profiles  which  are  in  good  agreement  with  observations 
and  numerical  model  results. 

The  solar  and  IR  parameterizations  are  summarized  in  the  following. 

a.  Short-Wave  Radiation 

Following  the  two-stream  model  solutions  by  Stephens  (Reference  42),  reflection  ( Re ), 
transmission  ( Tr )  and  absorption  (A)  are  expressed  as: 
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(25) 


Tr(fi0)  =  Au/R  (28) 

Ab(fJ,0)  =  1  —  ReifJ-o)  —  TrM  (29) 

where 

v?  =  [1  —  +  2/5(/i0)?50]/(  1  —  w0)  (^0) 

re//  “  {(1  ~  W0)[l  —  Wo  +  2/?(/i0)iyo]}^riV/ /^o  (^l) 

R  =  (u  +  l)2ea;p(re//)  -  (u  -  l)2exp(-re//)  (32) 

In  the  above  expressions,  fi0  is  the  zenith  angle,  r,v  is  the  optical  thickness  of  the  cloud, 
w0is  the  single-scattering  albedo,  and  &  is  the  backscattered  fraction  of  monodirectional  incident 
radiation.  The  values  for  w0  and  are  tabulated  as  functions  of  and  /i0- 

Now  solar  radiation  flux  profile  within  the  layered  cloud  is  computed  from  the  following 
equations  (Reference  39): 
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F(z)  =  FB  -  ( Fb  -  Fc){  1  -  exp{-{zB  -  z)/\s)}/ 
{1  -  exp[-(zB  -  2c)/As]} 

where 
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A 5  =  a(p0)W  +  b(p0)(  1  -  expi-^W  +  c(p0)}})  (34) 


Fc  =  Fb  +  AbFB  j 


(35) 


a(p0)  =  -0.022  +  0.038(1  -  p0) 


(36) 


b(p0)  =  56.8  -  14.7(1  -  p0)  (37) 

c(p0)  =  1.07-  1.15(1  -/zffl)  (38) 

In  the  above  expressions  W  is  the  total  cloud  liquid  water  content,  which  is  given  as, 

ZB 

W  =  j  p(z)Qc(z)dz  (39) 

where  p(z )  is  the  air  density,  zg  is  the  cloud  top,  and  zq  is  the  cloud  base.  The  units  of  W 
are  gm“2  and  7=0.021.  FB  =  -FB[{1-R(p0)}  where  FB[  is  the  downward  solar  flux  just  above 
the  top  of  the  cloud. 

b.  Long- Wave  Radiation 

This  parameterization  relies  on  values  of  external  conditions:  TB  and  Tc  are  the  absolute 
temperature  at  the  cloud  top  and  the  cloud  base,  respectively;  Gg|is  the  downward  IR  flux  at 
the  cloud  top;  Gclis  the  upward  IR  flux  at  the  cloud  base;  af=  0.13  and  a  [  =  0.158. 
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The  IR  flux  profiles  in  the  layered  cloud  are  given  as 

G{z)  —  Giexp[—(z  —  zc)I\l\  +  Guexp[—(zB  —  z)/\u\,  (40) 

where 


Gl  =  {^0  -  G\.exp[-(zB  -  zc)/\u]}/D, 

(41) 

Gu  =  {Gi  -  Goexp[-(zB  -  zc)/^l]}/D, 

(42) 

D  =  1  -  exp{-(zB  -  zc)/\n }, 

(43) 

1  _  i  1 

A n  A[/  \l 

(44) 

Go  =  Gc  T  -Bc  +  {Be  -  Gb  i)exp(-7]  |), 

(45) 

Gi  =  ( Gc  T  +Bb  -  2Bc)exp{-p  T)  +  Bb  -  GB  j 

(46) 

V  t,!=  ol  T,1  w, 

(47a, b) 

Be  =  <rT&, 

(48) 

Bb  =  aT4B, 

(49) 

\l  =  70W/(w~W*  +  2.GT), 

(50) 
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A v  =  140 W 


(51) 


—0.56 

where  a  is  the  Stefan-Boltzmann  constant. 

3.  Precipitation  Microphysics 

Following  Nickerson  (Reference  43),  the  microphysics  are  expressed  in  the  following 
equations: 


dQw  _  dQr \ 
ft  V  dt  )  sedi 


dQr  =  (dQA  (dQr\ 

V  dt  )  auto  V  dt  )  accr 


dQr)  {dQr) 

&  Jevap  V  dt  /  sedi 


(53) 


dNr  =  fdNr)  fdNr)  /  dNr) 

ft  V  ft  )  auto  V  dt  )  self  V  dt  )  sedi 


(54) 


where  Qw  is  the  sum  of  mixing  ratios  for  water  vapor,  cloud  water,  and  rain  water.  Qr  is 
the  rain  water  mixing  ratio  and  Nr  is  the  raindrop  number  concentration.  The  subscript  auto  is 
for  autoconversion  of  cloud  droplets  into  raindrops,  accr  is  for  accretion  of  cloud  droplets  by 
raindrops,  self  is  for  the  self-collection  of  the  raindrops,  evap  is  for  rain  evaporation  and  sedi  is 
for  rain  sedimentation.  Nickerson  (Reference  43)  discussed  in  detail  the  expressions  for  those 
microphysics  processes.  This  hydrological  model  is  substantially  more  complex  than  Kessler’s 
(Reference  44)  classical  parameterization  for  the  treatment  of  microphysics.  The  model  is  based 
on  a  more  realistic  log-normal  distribution  than  the  Marshall-Palmer  distribution  in  Kessler’s 
parameterization.  The  mean  diameter  of  raindrops  in  a  given  volume  is  estimated  from  equations 
for  the  mean  rain  water  mixing  ratio  and  the  raindrop  number  concentration.  The  terminal 
velocity  of  raindrops  in  each  volume  is  calculated  as  a  function  of  the  mean  raindrop  diameter, 
thus  the  terminal  velocity  values  in  the  model  are  expected  to  be  more  accurate  than  those 
obtained  by  a  simple  parameterization. 
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4.  Simulation-  of  Maritime  Stratus  Clouds 

A  numerical  simulation  was  conducted  to  investigate  complex  interactions  among  clouds, 
turbulence,  and  atmospheric  radiation.  Maritime  stratus  clouds  were  simulated  by  using  a  one¬ 
dimensional  version  of  HOTMAC  (Reference  20). 

Integration  initiated  at  0000  LT  on  day  200  and  continued  for  48  hours.  The  results  for  the 
second  24-hour  simulations  are  presented  here.  Figure  22  shows  a  time-height  variation  of  the 
mixing  ratio  of  cloud  liquid  water.  Clouds  were  thick  during  the  nocturnal  period  and  occupied 
the  layer  between  1200  meters  and  1600  meters  above  the  sea  surface.  As  the  short-wave  solai 
radiation  heating  increased,  clouds  dissipated  from  the  lower  part  and  the  clouds  became  as  thin 
as  100  meters  by  1800  LT.  As  the  solar  heating  subsided,  the  cloud  thickness  increased  again. 
Figure  23  shows  the  diurnal  variations  of  the  mixing  ratio  of  rain  water.  Rain  was  produced 
in  the  cloud  but  evaporated  before  it  reached  the  surface  except  during  a  short  period  between 
0600  and  0800. 

Figures  24  and  25  show,  respectively,  the  short-wave  radiation  heating  rate  (°C/day)  and 
the  long-wave  radiation  cooling  rate  (°C/day).  The  short-wave  radiation  heating  was  maximum 
at  the  cloud  top  and  decreased  exponentially  within  the  cloud  (Reference  39).  Short-wave 
radiation  heating  (Figure  24)  was  offset  by  long-wave  radiation  cooling  at  the  cloud  top  (Figure 
25).  Figure  25  also  indicates  long-wave  radiation  warming  at  the  cloud  base. 

No  attempt  was  made  to  compare  simulations  with  observations. 

A  one-dimensional  version  of  HOTMAC  was  used  to  simulate  time  evolution  of  fog 
over  horizontally  homogeneous  terrain.  The  results  were  compared  with  data  taken  at  Cabaw 
meteorological  tower  in  the  Netherlands  (Reference  31).  Wind,  temperature,  and  the  mixing 
ratio  of  water  vapor  were  initialized  with  measurements.  Simulations  were  initiated  at  midnight 
of  day  215  and  continued  for  24  hours. 

Fog  began  to  form  almost  immediately,  as  the  air  temperature  decreased  due  to  cooling  at 
the  ground  (Figure  26).  The  height  of  the  fog  increased  as  the  air  temperature  near  the  fog 
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Figure  22.  Diurnal  Variation  of  the  Cloud  Water  Mi? 

Cloud.  Contour  From  0  to  128  with  an  Ii 
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Figure  23.  Similar  to  Figure  22,  except  for  the  Rain  Water  Mixing  Ratio  (g/g).  Contour  from  0 
to  720  with  an  Interval  of  40.  Labels  Scaled  by  1.0e06. 
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Figure  24.  Similar  to  Figure  22,  except  for  the  Short-Wave  Radiation  Cooling  Rate  (°C/day). 
Contour  From  0  to  180  with  an  Interval  of  10. 
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5.  Simulation  of  Fog 


top  decreased  due  to  long-wave  radiation  cooling.  The  mixing  ratio  of  cloud  water  reached  a 
maximum  value  around  0300. 

Sunrise  was  shortly  after  0500,  but  fog  did  not  dissipate  until  0700  because  the  amount 
of  solar  energy  reached  the  ground  was  reduced  considerably,  due  to  reflection  and  absorption 
by  the  fog. 

Figure  27  shows  the  modeled  and  observed  temperature  and  mixing  ratio  of  water  vapor 
at  1.1  meters  agl.  The  mixing  ratio  at  the  ground  was  assumed  to  be  saturated  while  fog  existed 
and  a  constant  after  fog  dissipated.  Simulations  were  in  good  agreement  with  observations. 

E.  TASK  5:  Develop  an  innovative  method  to  treat  plumes  that  are  not  neutrally  buoyant. 

Highly  buoyant  plumes  modify  the  wind  and  turbulence  distributions  of  the  ambient  flow.  It 
is  almost  impossible  to  parameterize  or  express  such  modifications  without  deploying  a  dynamic 
plume  model.  A  physically  correct  way  to  treat  nonneutrally  buoyant  plumes  is  to  incorporate 
plume  dynamics  into  HOTMAC.  However,  the  dynamic  plume  model  requires  considerable 
computer  time  and  is  not  practical  for  emergency  modeling.  This  is  a  dilemma,  for  one  must 
choose  between  being  accurate  (but  impractical)  and  being  practical  (but  inaccurate).  It  is 
necessary  to  adopt  a  plume  parameterization  scheme  (instead  of  a  dynamic  plume  model)  to 
meet  the  time  constraints  of  an  emergency  response  modeling  system. 

We  have  tried  two  approaches:  (1)  modified  HOTMAC  to  include  plume  dynamics  (a 
fire  code)  and  (2)  modified  RAPTAD  only  to  incorporate  the  buoyant  plume  parameterization 
proposed  by  Van  Dop  (Reference  45).  The  fire  code  was  previously  run  on  a  Cray  supercomputer 
but  has  never  been  run  on  a  Sun  workstation.  The  fire  code  and  graphics  programs  were  modified 
to  run  on  a  workstation.  The  model  produced  the  results  identical  to  those  produced  on  the 
supercomputer.  But,  as  expected,  run  time  was  too  long  to  be  practical  for  emergency  modeling. 
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The  second  approach,  to  modify  RAPTAD  to  include  buoyant  plume  effects,  has  the  advantage 
of  being  faster.  The  disadvantage  is  that  it  does  not  consider  the  effects  of  a  buoyant  plume  on 
the  ambient  airflows.  In  other  words,  the  ambient  airflows  are  considered  to  be  unchanged  even 
if  a  buoyant  plume  modifies  both  wind  and  temperature  distributions.  Thus,  this  approach  is 
justifiable  only  when  ambient  air  modifications  are  confined  to  a  small  area  and  the  magnitude 


of  any  modification  is  small. 

The  vertical  velocity  of  a  buoyant  plume  is  computed  from  the  Langevin  equation  of  motion 
for  a  homogeneous  and  stationary  turbulent  flow.  The  temperature  of  a  buoyant  plume  is  also 
assumed  to  be  computed  by  the  Langevin  equation.  Following  Van  Dop  (Reference  45),  vertical 
velocity  (W),  buoyancy  (£),  and  plume  height  (Z)  are  computed  from  the  following  equations. 


where 


dW  w  ,  ^ 

dt  Tw+  ’ 

(55) 

-  =  -JL-n2w, 

dt  Tb 

(56) 

dt 

(57) 

B=g-{Q-  0a), 

(58) 

/y2  9  dea 

T  dz' 

(59) 

Tw  =  Tb  =  A(t  +  t0 ), 

(60) 

CO 

II 

(61) 

t0  =  1  sec. 

(62) 
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In  the  above  equations,  0  and  0a  are  potential  temperatures  of  plume  and  ambient  air, 
respectively. 

Test  simulations  were  conducted  under  neutral  conditions  (N2  =  0)  and  the  results  were 
compared  with  analytical  solutions.  Figure  28  shows  time  variations  of  B,  W,  and  Z  where 
initial  buoyancy  of  1  m/s2  and  initial  vertical  velocity  of  0.1  m/s  were  used.  Simulations  and 
analytical  solutions  are  in  good  agreement. 

Simulations  were  repeated  under  stable  (N2  >  0)  and  unstable  (N2  <  0)  atmospheric  conditions. 
Figures  29  and  30  show  time  variations  of  B,  W,  and  Z  when  N2  =  (0.408)2  and  N2  =  -(0.408)2, 
respectively.  Initial  buoyancy  of  1  m/s2  and  initial  vertical  velocity  of  0  m/s  were  used. 

Finally,  simulations  were  conducted  for  a  plume  whose  initial  density  is  heavier  than  the 
ambient  air.  Initial  buoyancy  of  -1  m/s2  and  initial  vertical  velocity  of  0  m/s  were  used.  Figures 
31  and  32  show  time  variations  of  B,  W,  and  Z  when  N2  =  (0.408)2  and  N2  =  -(0.408)2, 
respectively. 

All  simulations  were  conducted  with  an  integration  time  step  of  0.1  seconds.  Accuracy 
decreased  rapidly  with  increases  in  the  integration  time  step.  RAPTAD  uses  a  time  step  of  10 
seconds,  but  a  small  time  step  was  necessary  to  simulate  initial  plume  rise  (drop)  accurately. 
Therefore,  a  time  step  which  increased  linearly  with  time  was  used:  the  time  step  was  0. 1  seconds 
initially  and  became  10  seconds  at  100  seconds  after  the  release.  Figure  33  shows  time  variations 
of  B,  W,  and  Z  where  N2  =  -(0.408)2,  B0  =  -1,  and  a  variable  time  step  were  used.  The  results 
are  in  good  agreement  with  those  in  Figure  32,  where  a  constant  time  step  of  0. 1  was  used. 

F.  TASK  6:  Develop  a  method  to  predict  concentration  variances,  which  are  used  to  determine 
how  a  required  confidence  level  affects  the  extent  of  a  critical  area. 

Relatively  short  time-averaging  values  are  required  for  predicting  concentrations  of  toxic 
materials.  Such  values  normally  exhibit  great  variations  in  time  and  space.  Thus,  predictions  of 
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e  28.  Time  Variation  of  Buoyancy  (m/s2).  Vertical  Velocity  (m/s),  and  Plume  Height  (m) 
in  the  Neutral  Condition.  Initial  Buoyancy  of  1  m/s2  and  Initial  Vertical  Velocity 
of  0  m/s  were  used. 
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not  only  the  mean  but  also  the  variance  of  concentrations  are  essential  to  determine  the  extent 
of  a  critical  area  where  concentration  values  exceed  specified  limits.  The  size  of  a  critical  area 
becomes  larger  if  a  higher  confidence  level  is  required. 

Natural  variations  of  concentration  could  be  large  even  with  the  use  of  30-  to  60-minute 
averaging  times,  which  are  used  for  the  primary  toxics  at  VAFB.  The  initial  scope  of  this  task 
was  narrow,  focusing  on  the  specific  question  of  how  a  required  confidence  level  affects  plume 

footprint  size. 

RAPTAD  used  an  integration  time  step  of  10  seconds,  but  concentrations  were  sampled  every 
120  seconds  to  reduce  computational  time.  The  current  version  of  RAPTAD  allows  sampling 
at  up  to  50  stations. 

Figures  34  through  36  show  examples  of  time  variations  of  mean  concentration  and  ratios  of 
standard  deviation  to  mean  values  which  are  averaged  over  a  time  period  of  15  minutes.  These 
figures  indicate  that  the  standard  deviation  of  concentration  could  be  as  large  as  several  times 
the  mean  concentration,  particularly  near  the  plume  edge.  Although  the  ratios  were  at  their 
minimum  along  the  plume  axis  and  increased  rapidly  toward  the  plume  edge,  both  mean  values 
and  standard  deviations  were  at  their  maximum  along  the  plume  axis  and  decreased  toward 
the  plume  edge.  Mean  values  decreased  more  rapidly  than  standard  deviations,  thus  the  ratios 
increased  moving  from  the  axis  to  the  edge.  Only  a  limited  number  of  concentration  fluctuation 
measurements  are  available  for  comparison.  The  best  measurements  were  obtained  by  using 
wind  tunnels.  Preliminary  comparisons  indicate  that  our  simulations  are  qualitatively  in  good 
agreement  with  wind  tunnel  data.^  Both  standard  deviations,  and  ratios  of  standard  deviation  to 
mean,  showed  characteristics  similar  to  those  found  in  measurements  where  measured  mean  wind 
and  turbulence  were  used  in  RAPTAD  computation.  No  effort  was  made  to  compare  simulations 
with  atmospheric  data  because  we  are  not  aware  of  the  existence  of  any  such  data. 


Lee,  J.T.,  (Personal  communication),  1990. 
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Figure  34.  Examples  of  Time  Variations  of  Mean  Co 
Deviation  to  Mean  Concentration  at  Static 


G.  TASK  7:  Perform  model  verification  runs  with  a  nested  grid  version  of  HOTMAC  on  a 
workstation . 

A  nested  grid  version  of  HOTMAC  had  been  used  with  a  Cray  supercomputer,  but  it  had 
never  been  used  with  a  workstation  because  it  would  require  too  much  computer  time.  However, 
with  the  rapid  advancement  of  workstation  performance,  a  nested  grid  version  of  HOTMAC  can 
now  be  ran  on  a  workstation. 

Test  simulations  were  performed  to  demonstrate  systematically  how  nested  grids  improved 
the  simulation  of  wind  distributions.  Three  nested  grids  were  used  in  a  control  run.  Horizontal 
grid  spacings  (in  both  X  and  Y  directions)  were  500  meters,  2000  meters,  and  8000  meters  for 
the  innermost  grid  (Grid  3),  intermediate  grid  (Grid  2),  and  the  outer  grid  (Grid  1),  respectively. 
The  number  of  grids  used  in  the  horizontal  directions  were  32  X  32,  16  X  16,  and  8  X  8, 
respectively  for  Grid  3,  Grid  2,  and  Grid  1.  Each  grid  had  16  vertical  levels. 

Simulations  began  at  0600  LT  on  a  mid-August  day  (arbitrarily  chosen)  and  continued  for 
24  hours.  A  two-way  nesting  method  was  used:  the  outer  grid  provided  boundary  values  to  the 
inner  grid,  and  the  inner  grid  updated  outer  grid  values  with  new  values  at  grid  points  common 
to  both  grids. 

Figure  37  shows  1500  LT  wind  vector  distributions  at  6  meters  agl  at  1500  LT  where  grid 
spacings  in  the  horizontal  directions  are  8  kilometers.  Areas  enclosed  by  dashed  lines  indicate 
nested  grids.  Figures  38  and  39  show  wind  distributions  for  Grid  2  and  Grid  3.  As  expected, 
more  variations  in  space  were  simulated  when  gnd  spacing  decreased  to  2  kilometers  (Figure 
38)  and  500  meters  (Figure  39).  Wind  direction  was  northwesterly  and  wind  speed  was  3  m/s 
in  the  levels  above  the  boundary  layer.  Wind  directions  in  the  surface  layer  show  considerable 
variations  due  to  nonhomogeneous  temperature  distributions  over  complex  terrain.  Sea-breeze 
and  upslope  flows  were  clearly  seen  in  Grid  3  (Figure  39). 

Wind  directions  in  the  surface  layer  vary  considerably  with  time.  Figure  40  shows  next 
morning  0200  LT  wind  distributions  at  6  meters  agl  in  the  Grid  1  computational  domain  (64 
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Figure  37.  Wind  Distributions  at  6  Meters  agl  at  1500  LT  on  Day  233.  Arrows  indicate  Wind 
Vectors,  whose  Scale  is  shown  at  the  Left  Upper  Comer.  Solid  Lines  indicate  Ground 
Elevation  Contours  in  increments  of  50  Meters.  Areas  enclosed  by  Dashed  Lines 
show  Nested  Grids.  Computation  Domains  are  referred  to  as  Grid  1,  Gnd  2,  and 
Grid  3  for  the  Outer,  Intermediate,  and  Inner  Grids,  respectively. 
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kilometers  X  64  kilometers),  where  the  grid  spacing  was  8  kilometers.  Figures  41  and  42  show 
corresponding  wind  distributions  over  the  Grid  2  and  Grid  3  computational  domains.  Sea  breezes 
were  replaced  by  land  breezes,  and  upslope  flows  were  replaced  by  downslope  flows.  These 
variations  are  common  in  coastal  complex-terrain  areas,  provided  that  synoptic  wind  forcing 
is  small. 

Figures  43  through  46  show  wind  distributions  where  only  two  grids  were  used,  in  comparison 
with  the  control  run,  where  three  grids  were  used.  The  innermost  grid  of  the  control  run  was 
omitted.  Figures  43  and  44  are  control  run  counterparts  of  Figures  37  and  38,  respectively. 
Atmospheric  turbulence  was  significant,  due  to  strong  heating  from  the  ground  and  enhanced 
momentum  transfer  in  the  vertical  direction.  This  resulted  in  relatively  uniform  wind  distributions 
in  space.  On  the  other  hand,  turbulence  mixing  was  at  a  minimum  during  the  nocturnal  period, 
due  to  stable  density  stratification  created  by  radiation  cooling  at  the  ground.  Thus,  winds  in 
the  surface  layer  were  often  decoupled  from  the  winds  in  the  upper  levels.  Figures  45  and  46 
are  0200  LT  wind  distributions  at  6  meters  agl  where  two  grids  were  used.  Figures  45  and  46 
are  control  run  counterparts  of  Figures  40  and  41.  There  are  noticeable  differences  in  the  wind 
distributions  between  the  three-grid  (control)  and  two-grid  runs,  particularly  in  the  areas  where 
Grid  3  was  nested  in  the  control  run. 

Finally,  simulations  were  conducted  without  a  nested  grid.  Figures  47  and  48  show  1500  LT 
and  0200  LT  wind  distributions  at  6  meters  agl.  Wind  distributions  at  1500  LT  (Figure  47)  are 
similar  to  those  for  the  control  run  (Figure  37)  because  of  strong  mixing  in  the  vertical  direction, 
due  to  turbulence,  which  resulted  in  almost  uniform  wind  distributions.  On  the  other  hand,  wind 
distributions  at  0200  LT  (Figure  48)  are  considerably  different  from  those  for  the  control  run 
(Figure  40),  particularly  in  the  area  where  Grid  3  was  nested. 
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Figure  42.  Similar  to  Figure  40  except  showing  Grid  3. 
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Figure  44.  Similar  to  Figure  38,  but  using  only  Grid  1  and  Grid  2 
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Figure  46.  Similar  to  Figure  41,  but  using  only  Grid  1  and  Grid  2 


H.  TASK  8:  Upgrade  computer  capabilities  used  for  emergency  response  modeling  at  VAFB. 


HOTMAC  and  RAPTAD  have  been  run  on  several  workstations  produced  by  leading  computer 
companies,  including  Sun  Microsystems,  Silicon  Graphics,  Data  General,  Hewlett-Packard,  and 
Digital  Data  Corporation.  We  have  developed  a  user-friendly  menu  system  which  includes  a 
variety  of  graphic  displays  for  the  Sun  Microsystems  and  Silicon  Graphics  workstations.  We 
rewrote  our  interface  and  graphics  software  for  the  Sun’s  OpenLook  window. 

We  performed  the  first  half  of  this  project  with  YSA  computer  equipment  in  order  to:  (1) 
fully  coordinate  the  feasibility  and  support  for  installation  of  a  suggested  workstation  at  VAFB, 
and  (2)  take  advantage  of  any  performance  gains  which  might  become  available  during  the  first 
half  of  the  project  period. 

Based  on  discussions  between  the  personnel  at  VAFB  and  YSA,  the  following  computer 
hardware  and  software  were  selected: 

1.  A  Sun  SPARCstation  10  with  32  Mb  RAM,  19"  color  monitor,  400  Mb  internal  SCSI 
disk,  1.44  Mb  3l”  internal  floppy  disk. 

2.  A  1.2  Gb  external  disk. 

3.  A  150  Mb  tape  drive. 

4.  A  CD  drive. 

5.  Sun  FORTRAN  compiler. 

6.  Sun  C  compiler. 

7.  NCAR  Graphics. 

Hardware  and  software  were  installed  and  tested  at  YSA,  and  delivered  to  VAFB. 
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SECTION  IV 


CONCLUSION  AND  RECOMMENDATIONS 


A.  CONCLUSION 

•  It  is  feasible  to  operate,  on  a  workstation,  a  nested-grid,  three-dimensional  atmospheric 
model,  HOTMAC,  and  a  three-dimensional  dispersion  model,  RAPTAD,  to  forecast  the 
transport  and  diffusion  of  airborne  materials  at  VAFB. 

•  The  computer  capabilities  for  emergency  response  applications  at  VAFB  were  upgraded. 
HOTMAC  and  RAPTAD  were  installed  and  tested  on  a  Sun  SPARCstationlO. 

•  HOTMAC  was  modified  to  run  continuously  and  store  wind  and  turbulence  forecasts  for 
the  next  24  hours.  When  an  emergency  occurs,  RAPTAD  can  be  used  immediately  with 
the  wind  data  stored  on  disk.  The  RAPTAD  computation  is  much  faster  than  that  of 
HOTMAC.  Thus,  this  approach  meets  better  the  time  constraints  of  emergency  situations. 

•  A  method  was  developed  to  integrate  large-scale  weather  data  into  HOTMAC.  The 
weather  data  are  used  to  initialize  and  correct  HOTMAC  forecasts.  A  four-dimensional 
data  assimilation  method  was  used. 

•  A  method  was  developed  to  predict  concentration  variances  which  can  be  used  to  estimate 
uncertainties  associated  with  predictions. 

•  Model  physics  was  added  to  simulate  fog  formation  and  dissipation.  Fog  is  frequently 
observed  at  VAFB  and  affects  the  heat  energy  balance  at  ground  level. 

•  Positive  and  negative  buoyancy  effects  of  plumes  were  incorporated  into  RAPTAD. 

B.  RECOMMENDATIONS 

•  Become  familiar  with  the  background  theories  and  operating  procedures  of  HOTMAC 
and  RAPTAD.  Increase  experience  by  running  the  models  under  different  weather 
conditions. 
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•  Consider  customizing  pre-  and  post-  processors  of  HOTMAC  and  RAPTAD  to  meet 
specific  needs  at  VAFB.  The  current  versions  of  HOTMAC  and  RAPTAD  are  designed 
for  general  applications. 

•  Continue  to  upgrade  computer  capabilities  at  VAFB.  Computer  hardware  technology  is 
expected  to  advance  further  and  to  allow  much  more  sophisticated  models  to  be  operated 
in  much  less  time  than  previously  considered  possible. 

•  Continue  to  improve  the  model  physics  of  HOTMAC  and  RAPTAD.  These  codes  have 
model  physics  which  is  considered  to  be  the  state-of-the-science,  but  falls  short  of 
replicating  the  complex  physics  in  the  real  atmosphere.  In  particular,  we  strongly 
recommend  that  the  new  model  physics  added  here  (precipitation  microphysics,  radiation 
transfer  in  fog  and  clouds,  and  positive  and  negative  plume  buoyancy  effects)  be  tested 
further  by  comparing  simulations  with  observations. 
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